reconstruct.cpp

Go to the documentation of this file.

/* reconstruct.cpp : Sample reconstruction program
 * Copyright (C) 2015  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 <iostream>

#include "atomprobe/atomprobe.h"

using namespace std;
using namespace AtomProbe;


int main( int argc, char *argv[])
{
    string filename;    
    if(argc ==1)
    {
        cout << "Enter file name" << endl;

        cin >> filename;
    }
    else if(argc ==2)
        filename=argv[1];
    else
    {
        cerr<< "USAGE: " << argv[0] << " [FILENAME]" << endl;
        return 1;
    }

    vector<vector<float> > data;
    vector<string> headers;
    cerr << "Loading :" << filename << endl;
    if(loadTextData(filename.c_str(),data,headers,"\t",true,true,1))
    {
        cerr << " Failed to load text data" << endl;
        return 1;
    }

    //Find the header that is det_x, det_y
    vector<string> search = { "det_x","det_y"};
    vector<unsigned int> offset;
    offset.resize(search.size(),(unsigned int)-1);
    for(auto ui=0u; ui<headers.size();ui++)
    {
        auto s = std::find(search.begin(),search.end(),
                stripWhite(headers[ui]));

        if(s != search.end())
            offset[distance(search.begin(),s)] = ui;
    }

    for(auto h : headers)
        cerr << "Header :" << h << endl;

    if(find(offset.begin(),offset.end(),-1) != offset.end())
    {
        cerr << "Unable to locate det_x or det_y field" << endl;
        return 1;
    }


    vector<float> detX,detY;
    detX.swap(data[offset[0]]);
    detY.swap(data[offset[1]]);
    data.clear();

    cerr << "Loaded :" << detX.size() << " lines." << endl;
    cerr << "Loaded :" << detY.size() << " lines." << endl;
    float detBounds[2][2];
    detBounds[0][0]=*std::min_element(detX.begin(),detX.end());
    detBounds[1][0]=*std::min_element(detY.begin(),detY.end());
    detBounds[0][1]=*std::max_element(detX.begin(),detX.end());
    detBounds[1][1]=*std::max_element(detY.begin(),detY.end());

    cerr << "Pre-cut" <<endl;
    cerr <<  detBounds[0][0]<< "<->" << detBounds[0][1] << endl;
    cerr <<  detBounds[1][0]<< "<->" << detBounds[1][1] << endl;
    


    const float DET_RADIUS=0.04; //meters
    cerr << "Trimming to radius :" << DET_RADIUS << endl;
    vector<float> trimX,trimY;
    for(auto ui=0u;ui<detX.size();ui++)
    {
        if(detX[ui]*detX[ui] + detY[ui]*detY[ui] < DET_RADIUS*DET_RADIUS)
        {
            trimX.push_back(detX[ui]);
            trimY.push_back(detY[ui]);
        }
    }

    detX.swap(trimX);
    detY.swap(trimY);

    AtomProbe::ReconstructionSphereOnCone r;
    AtomProbe::ModifiedFocusSphericProjection model(1.05f);

    std::vector<float> ionV,tof; //Ion volumes

    //Create initial parameters
    r.setInitialRadius(35);
    r.setProjModel(&model); //Use the Modified Stereographic model
    r.setReconFOV(atan(DET_RADIUS/0.12));  //Set the 1/2 FOV
    r.setFlightPath(0.12); //Set flight path
    r.setShankAngle(15.2*M_PI/180.0f); //Angle as measured from recon image
    r.setDetectorEfficiency(1); //Set efficiency

    //make space for ion volumes & TOF
    ionV.resize(detX.size()); 
    tof.resize(detX.size());

    //OK, so now we know which column we want
    using AtomProbe::IonHit;

    std::vector<IonHit> pts;

    bool errorRes;
    errorRes = r.reconstruct(detX,detY,tof,ionV,pts);

    if(!errorRes)
    {
        cerr << "Reconstruction failed for reasons" << endl;
    }
    else
    {
        cerr << "Successfully reconstructed!" << endl;
        if(savePosFile(pts,"reconstructed.pos"))
        {
            cerr << "Error writing reconstruction to file" << endl;
        }
    }
}