pulse-filter.cpp

Go to the documentation of this file.

/* pulse-filter.cpp : Implements the pulse filtering method of Yao et al, MethodsX, 2016; 3: 268–273.  
 * Copyright (C) 2017  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include <atomprobe/atomprobe.h>

#include <gsl/gsl_histogram2d.h>


#include <iostream>
#include <cstdlib>
#include <vector>
#include <fstream>


using std::cerr;
using std::endl;
using std::vector;
using namespace AtomProbe;

enum
{
    ERR_BAD_ARGS,
    ERR_BAD_EPOS_READ,
    ERR_BAD_NUMSAMPLES,
    ERR_NOMEM
};


bool convertToPos(const vector<EPOS_ENTRY> &epos, vector<IonHit> &pos)
{
    try
    {
        pos.resize(epos.size());
    }
    catch(std::bad_alloc)
    {
        return false;
    }
#pragma omp parallel for
    for(size_t ui=0;ui<epos.size();ui++)
    {
        pos[ui].setPos(epos[ui].xDetector,
                 epos[ui].yDetector,
                    epos[ui].z);
        pos[ui].setMassToCharge( epos[ui].hitMultiplicity);
    }

    return true;
}

bool dumpHistogram(gsl_histogram2d *h, const char *outFile)
{
    std::ofstream of(outFile);

    if(!of)
        return false;

    for(unsigned int ui=0;ui<h->nx; ui++)
    {
        for(unsigned int uj=0;uj<h->ny; uj++)
        {
            of<< gsl_histogram2d_get(h,ui,uj) << " ";
        }
        of << endl;
    }

    return true;
}


void filterEposByPulse(const vector<EPOS_ENTRY> &input, unsigned int filterCount,
                        vector<EPOS_ENTRY> &output)
{
    for(unsigned int ui=0;ui<input.size(); ui++)
    {
        //If we reach an ion with a small deltaPulse count,
        // then we retain it. a deltaPulse of zero means its all the  
        // same ion
        while(ui < input.size() &&
            input[ui].deltaPulse <filterCount) 
        {
            output.push_back(input[ui]);
            ui++;
        }
    }
}


int main(int argc, const char *argv[])
{

    if(argc !=4)
    {
        cerr << "USAGE: " << argv[0] << " EPOSFILE SAMPLES OUTPUTPOS" << endl;
        cerr << endl;
        cerr << " This filters an epos file by the number of pulses between evaporation events. Only events with less than SAMPLES pulses between them are kept" << endl;

        cerr << "See Yao et al, MethodsX, 2016, 3:268-273" << endl;
        return ERR_BAD_ARGS;
    }


    unsigned int filterCount;   
    if(stream_cast(filterCount,argv[2]))
    {
        cerr << "Unable to interpret filter count value" << endl;
        return ERR_BAD_ARGS;
    
    }

    if(!filterCount)
    {
        cerr << "filter count must be positive" << endl;
        return ERR_BAD_ARGS;
    }

    cerr << "Filtering EPOS file....";
    unsigned int CHUNK_SIZE=4096;

    ProgressBar pb;
    pb.setLength(50);
    pb.init();

    unsigned int outputSize=0,chunkOffset=0,totalEntries=0,entriesLeft;
    vector<EPOS_ENTRY> eposEntries,eposOutput;
    EPOS_ENTRY lastEntry;
    do
    {
        vector<EPOS_ENTRY> thisEposOutput;
    
    
        unsigned int errCode;
        if((errCode=chunkLoadEposFile(eposEntries,argv[1],
            CHUNK_SIZE,chunkOffset,entriesLeft)))
        {
            cerr << endl << "\tError reading epos, " << argv[1]  << " : " << RECORDREAD_ERR_STRINGS[errCode] << endl;
            return ERR_BAD_EPOS_READ;
        }

        //include the tail entry from last run, if applicable
        if(chunkOffset)
        {
            //FIXME:Unclear if this is correct. We would need
            // all of the entries which were not nMultiples==0? 
            eposEntries.insert(eposEntries.begin(),lastEntry);
        }

        //Update our offset
        chunkOffset++;

        //Record the last entry, as we need it later
        lastEntry = eposEntries.back();

        //Filter the epos files 
        totalEntries+=eposEntries.size();
        filterEposByPulse(eposEntries,filterCount,thisEposOutput);

        //If there is more to process, do not include the 
        // last entry in the processing chain in our results, as we need 
        // it to fix the "seam" in the chain of pulses of epos 
        // files, by appending it before filtering

        if(thisEposOutput.size() && 
            (thisEposOutput.back() == lastEntry) && entriesLeft)
            thisEposOutput.pop_back();

        //Convert epos to pos
        vector<IonHit> ionHits;
        ionHits.resize(thisEposOutput.size());

        for(auto ui=0;ui<thisEposOutput.size();ui++)
            thisEposOutput[ui].getIonHit(ionHits[ui]);

        //Append to the pos file
        const bool APPEND=true;
        if(savePosFile(ionHits,argv[3],APPEND))
        {
            cerr << "Error writing pos file" << endl;
            return 2;

        }

        outputSize+=thisEposOutput.size();
        thisEposOutput.clear();


        //Update progress
        pb.update(100.0f*(float)chunkOffset/(float)(chunkOffset+entriesLeft/CHUNK_SIZE));
    }while(entriesLeft);

    //Finalise progress bar
    pb.finish();

    //Report back to user
    cerr << "\tFiltered " << totalEntries << " entries." << endl;
    cerr << outputSize<< " entries after filtering." << endl;

    

}