reconstruction-simple.cpp

Go to the documentation of this file.

/* reconstruction-simple.cpp :  Mass overlap deconvolution algorithm
 * 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/reconstruction/reconstruction-simple.h"

#include "atomprobe/helper/aptAssert.h"
#include "atomprobe/primitives/boundcube.h"
#include "atomprobe/io/dataFiles.h"

#include "atomprobe/helper/constants.h"

#include <vector>
#include <random>

namespace AtomProbe{

ReconstructionSphereOnCone::ReconstructionSphereOnCone()
{ 
    reconFOV=initialRadius=shankAngle=0;
    detEfficiency=1;
    evolutionMode=EVOLUTION_MODE_SHANKANGLE;
}

void ReconstructionSphereOnCone::setProjModel(SphericPlaneProjection *model)
{
    projectionModel=model;
}

void ReconstructionSphereOnCone::setRadiusEvolutionMode(unsigned int mode)
{
    evolutionMode=mode;
}

void ReconstructionSphereOnCone::setDetectorEfficiency(float detEff)
{
    ASSERT(detEfficiency > 0 && detEfficiency <=1);
    detEfficiency =detEff;
}

void ReconstructionSphereOnCone::setFlightPath(float fp)
{
    ASSERT(fp> 0);
    flightPath=fp;
}

void ReconstructionSphereOnCone::setReconFOV(float angle)
{
    ASSERT(angle< M_PI)
    reconFOV=angle;
}

void ReconstructionSphereOnCone::setShankAngle(float angle)
{
    ASSERT(shankAngle >=0 && shankAngle < M_PI/2.0f);
    shankAngle=angle;
}


void ReconstructionSphereOnCone::setInitialRadius(float r)
{
    ASSERT(initialRadius >=0);
    initialRadius=r;
}

bool ReconstructionSphereOnCone::reconstruct(const std::vector<float> &detectorX, 
                        const std::vector<float> &detectorY,
                        const std::vector<float> &tof,
                        const std::vector<float> &ionVolume,
                        std::vector<IonHit> &outputPts) const
{
    ASSERT(evolutionMode == EVOLUTION_MODE_SHANKANGLE);
    ASSERT(detectorX.size() == detectorY.size());
    ASSERT(detectorX.size() == ionVolume.size());   
    outputPts.resize(detectorX.size());

    
    //Distance from sphere apex to cone apex
    double initZoff;

    //This is for the tangent to the sphere taking angle shankangle
    //===========

    if(shankAngle > sqrt(std::numeric_limits<float>::epsilon()))
        initZoff= initialRadius*(1.0/sin(shankAngle) -1.0) ;
    else
        initZoff=initialRadius;
    //============
   
    float etaFOV =  projectionModel->thetaToEta(reconFOV);

    //This is the initial height from origin to edge of cone
    double distRatio=initialRadius/initZoff;
    

    //set current Z offset to initial Z offset to sphere
    double zOffset=initZoff;
    //Set current radius
    double currentR;


    //Reconstruction is done in two passes to speed it up.
    // the first pass cannot be readily parallelised, but the second can.
    
    //=== PASS 1 ==
    //Switch between dynamic and static radius loops
    // First just compute the Z position that corresponds to the centre of the recon sphere
    // we do this so we can separate non-parallel and parallel regions
    if(shankAngle > std::numeric_limits<float>::epsilon())
    {
        //This, divided by radius^2, gives inverse area of spherical cap
        float areaFactor = 1.0f/(2.0f*M_PI*(1.0-cos(reconFOV)));
    
        for(unsigned int ui=0;ui<detectorX.size() ;ui++)
        {

            //Use similar triangles to scale-up size
            //R'  = R/L*(L +x), where (L+x) is the current Z offset, 
            //  L is the initial offset, and R is initial radius
            // See doc/figures/shank-reconstruction.svg for diagram
            currentR= distRatio*(zOffset);

            //set the delta from the cone apex
            outputPts[ui][2]=zOffset;
        
            //Update Z-offset by shuffling by spherical cap area    
            zOffset +=(ionVolume[ui]/detEfficiency)*areaFactor/(currentR*currentR);

        }
    }
    else
    {
        //We can compute this exactly once, as it wont change
        float areaSphericalCapInv;
        areaSphericalCapInv = 1.0/(2.0*M_PI*initialRadius*initialRadius*(1.0-cos(reconFOV)));

        
        //TODO: Parallel reduction sum
        for(unsigned int ui=0;ui<detectorX.size() ;ui++)
        {
            //set the origin of the recon sphere, which we will later update with complete point
            outputPts[ui][2]=(zOffset-initZoff);
        
            //Update Z-offset
            zOffset +=(ionVolume[ui]/detEfficiency)*areaSphericalCapInv;

        }
    }
    //=============

    //=== PASS 2===
    //Expand point to the full position on sphere
#ifdef _OPENMP
    bool spin=false;
#endif

    #pragma omp parallel for
    for(unsigned int ui=0;ui<detectorX.size() ;ui++)
    {
#ifdef _OPENMP
        if(spin)
            continue;
#endif

        float px,py ;
        projectionModel->scaleDown(flightPath,detectorX[ui],
                detectorY[ui],px,py);   

        float eta,phi;
        if(!projectionModel->toAzimuthal(px,py,eta,phi))
        {
#ifndef _OPENMP
            return false;
#else
            spin=true;
#endif
        }

        currentR= distRatio*(outputPts[ui][2]);
        outputPts[ui][0]=currentR*sin(eta)*cos(phi);
        outputPts[ui][1]=currentR*sin(eta)*sin(phi);
        outputPts[ui][2]+=currentR*(1.0f-cos(eta)) - initZoff;


        outputPts[ui][3]=tof[ui];
    }

#ifdef _OPENMP
    if(spin)
        return false;
#endif
    //=============

    return true;
}


#ifdef DEBUG

bool reconstructTest()
{
    //The goal of this is to "exercise" the reconstruction code
    //Create a random "detector sequence"

    std::vector<float> dx, dy, tof,ionV;

    const unsigned int NUM_PTS=100000;

    //Ion volume is very exaggerated
    ionV.resize(NUM_PTS,0.1);
    dx.resize(NUM_PTS);
    dy.resize(NUM_PTS);
    tof.resize(NUM_PTS);

    std::random_device r;
    std::mt19937 mtRnd(r());
    std::uniform_real_distribution<float> ud(0,1.0f);

    //detector total width
    const float DET_SIZE=40.0f; //mm
    const float FLIGHT_PATH=42.0f; //mm
    for(unsigned int ui=0;ui<NUM_PTS;ui++)
    {
        do
        {
            dx[ui] = 2.0*(ud(mtRnd)-0.5f); //Normalised -1->1 space
            dy[ui] = 2.0*(ud(mtRnd)-0.5f);
        } while(dx[ui]*dx[ui] + dy[ui]*dy[ui] > 1.0f);

        //Convert to detector size
        // Divide by two,as the detector spans 0-1 sapace
        dx[ui]*=DET_SIZE/2.0f;
        dy[ui]*=DET_SIZE/2.0f;

        tof[ui]=ud(mtRnd)+0.5;
    }



    ReconstructionSphereOnCone recons;

    //Specify the model for reconstruction. 
    // Gnomonic is a sphere-centred projection  
    StereographicProjection model;
    recons.setProjModel(&model);

    const float TIP_RADIUS=15.0f;
    //Provide reconstruction parameters
    recons.setInitialRadius(TIP_RADIUS);
    recons.setReconFOV(atan(0.5*DET_SIZE/FLIGHT_PATH));
    recons.setShankAngle(10.0f*M_PI/180.0f);
    recons.setFlightPath(FLIGHT_PATH);
    std::vector<IonHit> pts;
    bool errRes;
    errRes=recons.reconstruct(dx,dy,tof,ionV,pts);

    TEST(errRes,"Reconstruction should not fail");
    TEST(pts.size() && (pts.size() <= NUM_PTS),"One event per reconstruction");
    
    BoundCube bc;
    bc.setBounds(pts);
    TEST(bc.isValid(),"Valid bounding cube");
    //Boundcube should have a nonzero volume
    TEST(bc.getVolume() > 1e-5,"Recon should have some volume");

    const double FUDGE=1.5;
    TEST(bc.getSize(0) <= 2.0*TIP_RADIUS*FUDGE,
        "Bounding box should not exceed final tip diameter");

    TEST(bc.getSize(1) <= 2.0*TIP_RADIUS*FUDGE,
        "Bounding box should not exceed final tip diameter");
    
    savePosFile(pts,"reconstruct-test.pos");
    return true;    
}
#endif
}