cpu_tracker.h

Go to the documentation of this file.

#pragma once

#include "QueuedTracker.h"
#include "utils.h"
#include "scalar_types.h"
#include "kissfft.h"

class XCor1DBuffer {
public:
    XCor1DBuffer(int xcorw) : fft_forward(xcorw, false), fft_backward(xcorw, true), xcorw(xcorw)
    {}

    kissfft<scalar_t> fft_forward;  
    kissfft<scalar_t> fft_backward; 
    int xcorw;                      

    void XCorFFTHelper(complex_t* prof, complex_t* prof_rev, scalar_t* result);
};

class CPUTracker
{
public:
    int width;          
    int height;         
    int xcorw;          
    int trackerID;      

    float *srcImage;    
    float *debugImage;  
    float mean;         
    float stdev;        
#ifdef _DEBUG
    float maxImageValue;    
#endif

    float* zluts;
    bool zlut_memoryOwner;              
    int zlut_planes;                    
    int zlut_res;                       
    int zlut_count;                     
    float zlut_minradius;               
    float zlut_maxradius;               
    std::vector<float> zlut_radialweights;  
    kissfft<scalar_t> *qa_fft_forward;  
    kissfft<scalar_t> *qa_fft_backward; 

    bool testRun;

    float* GetRadialZLUT(int index)  { return &zluts[zlut_res*zlut_planes*index]; }

    XCor1DBuffer* xcorBuffer;                   
    std::vector<vector2f> quadrantDirs;         
    int qi_radialsteps;                         
    kissfft<scalar_t> *qi_fft_forward;          
    kissfft<scalar_t> *qi_fft_backward;         

    class FFT2D {
    public:
        kissfft<float> xfft;                    
        kissfft<float> yfft;                    
        std::complex<float> *cbuf;              

        FFT2D(int w, int h) : xfft(w,false), yfft(h,false) {cbuf=new std::complex<float>[w*h]; }
        
        ~FFT2D() { delete[] cbuf; }

        void Apply(float* d);
    };
    FFT2D *fft2d;               

    float& GetPixel(int x, int y) { return srcImage[width*y+x]; }   
    int GetWidth() { return width; }                                
    int GetHeight() { return height; }                              

    CPUTracker(int w, int h, int xcorwindow=128, bool testRun = false);

    ~CPUTracker();

    bool CheckBoundaries(vector2f center, float radius);

    vector2f ComputeXCorInterpolated(vector2f initial, int iterations, int profileWidth, bool& boundaryHit);

    vector2f ComputeQI(vector2f initial, int iterations, int radialSteps, int angularStepsPerQuadrant, 
        float angStepIterationFactor, float minRadius, float maxRadius, bool& boundaryHit, float* radialweights=0);

    struct Gauss2DResult {
        vector2f pos;           
        float I0;               
        float bg;               
    };

    Gauss2DResult Compute2DGaussianMLE(vector2f initial, int iterations, float sigma);

    scalar_t QI_ComputeOffset(complex_t* qi_profile, int nr, int axisForDebug);

    scalar_t QuadrantAlign_ComputeOffset(complex_t* profile, complex_t* zlut_prof_fft, int nr, int axisForDebug);

    float ComputeAsymmetry(vector2f center, int radialSteps, int angularSteps, float minRadius, float maxRadius, float *dstAngProf=0);

    template<typename TPixel> void SetImage(TPixel* srcImage, uint srcpitch);

    void SetImage16Bit(ushort* srcImage, uint srcpitch) { SetImage(srcImage, srcpitch); }

    void SetImage8Bit(uchar* srcImage, uint srcpitch) { SetImage(srcImage, srcpitch); }

    void SetImageFloat(float* srcImage);

    void SaveImage(const char *filename);

    vector2f ComputeMeanAndCOM(float bgcorrection=0.0f);

    void ComputeRadialProfile(float* dst, int radialSteps, int angularSteps, float minradius, float maxradius, vector2f center, bool crp, bool* boundaryHit=0, bool normalize=true);
    
    void ComputeQuadrantProfile(scalar_t* dst, int radialSteps, int angularSteps, int quadrant, float minRadius, float maxRadius, vector2f center, float* radialWeights=0);
    
    float ComputeZ(vector2f center, int angularSteps, int zlutIndex, bool* boundaryHit=0, float* profile=0, float* cmpprof=0, bool normalizeProfile=true)
    {
        float* prof = profile ? profile : ALLOCA_ARRAY(float, zlut_res);
        ComputeRadialProfile(prof,zlut_res,angularSteps, zlut_minradius, zlut_maxradius, center, false, boundaryHit, normalizeProfile);
        return LUTProfileCompare(prof, zlutIndex, cmpprof, LUTProfMaxQuadraticFit);
    }
    
    void FourierTransform2D();

    void FourierRadialProfile(float* dst, int radialSteps, int angularSteps, float minradius, float maxradius);

    void Normalize(float *image=0);

    void SetRadialZLUT(float* data, int planes, int res, int num_zluts, float minradius, float maxradius, bool copyMemory, bool useCorrelation);
    
    void SetRadialWeights(float* radweights);

    enum LUTProfileMaxComputeMode { 
        LUTProfMaxQuadraticFit,         
        LUTProfMaxSplineFit,            
        LUTProfMaxSimpleInterp          
    };

    void CalculateErrorCurve(double* errorcurve_dest, float* profile, float* zlut_sel);

    void CalculateInterpolatedZLUTProfile(float* profile_dest, float z, int zlutIndex);

    float CalculateErrorFlag(double* prof1, double* prof2);

    float LUTProfileCompare(float* profile, int zlutIndex, float* cmpProf, LUTProfileMaxComputeMode maxPosMethod, float* fitcurve=0, int *maxPos=0, int frameNum = 0);
    
    float LUTProfileCompareAdjustedWeights(float* rprof, int zlutIndex, float z_estim);

    float* GetDebugImage() { return debugImage; }

    void ApplyOffsetGain(float *offset, float *gain, float offsetFactor, float gainFactor);

    void AllocateQIFFTs(int nsteps);

    vector3f QuadrantAlign(vector3f initial, int beadIndex, int angularStepsPerQuadrant, bool& boundaryHit);
};

template<typename TPixel>
void CPUTracker::SetImage(TPixel* data, uint pitchInBytes)
{
    uchar* bp = (uchar*)data;

    for (int y=0;y<height;y++) {
        for (int x=0;x<width;x++) {
            srcImage[y*width+x] = ((TPixel*)bp)[x];
        }
        bp += pitchInBytes;
    }

    mean=0.0f;
}