CPUTracker Class Reference

Class with all CPU algorithm implementations. More...

#include <cpu_tracker.h>

Classes
class	FFT2D
	Class to facilitate 2D Fourier transforms. More...
struct	Gauss2DResult
	Structure to group results from the 2D Gaussian fit. More...

Public Types
enum	LUTProfileMaxComputeMode { LUTProfMaxQuadraticFit, LUTProfMaxSplineFit, LUTProfMaxSimpleInterp }
	Settings to compute the Z coordinate from the error curve. More...

Public Member Functions
float *	GetRadialZLUT (int index)
	Get the start of the ZLUT of a specific bead. More...
float &	GetPixel (int x, int y)
	Get the image pixel greyscale value at point (x,y). More...
int	GetWidth ()
	Get the width of the image. More...
int	GetHeight ()
	Get the height of the image. More...
	CPUTracker (int w, int h, int xcorwindow=128, bool testRun=false)
	Create an instance of CPUTracker. More...
	~CPUTracker ()
	Destroy a CPUTracker instance. More...
bool	CheckBoundaries (vector2f center, float radius)
	Test to see if a circle extends outside of image boundaries. More...
vector2f	ComputeXCorInterpolated (vector2f initial, int iterations, int profileWidth, bool &boundaryHit)
	Compute the cross correlation offsets and resulting position. More...
vector2f	ComputeQI (vector2f initial, int iterations, int radialSteps, int angularStepsPerQuadrant, float angStepIterationFactor, float minRadius, float maxRadius, bool &boundaryHit, float *radialweights=0)
	Execute the quadrant interpolation algorithm. More...
Gauss2DResult	Compute2DGaussianMLE (vector2f initial, int iterations, float sigma)
	Calculate a 2D Gaussian fit on the image. More...
scalar_t	QI_ComputeOffset (complex_t *qi_profile, int nr, int axisForDebug)
	QI helper function to calculate the offset from concatenated profiles. More...
scalar_t	QuadrantAlign_ComputeOffset (complex_t *profile, complex_t *zlut_prof_fft, int nr, int axisForDebug)
	QA helper function to calculate the offset from concatenated profiles. More...
float	ComputeAsymmetry (vector2f center, int radialSteps, int angularSteps, float minRadius, float maxRadius, float *dstAngProf=0)
	Find a measure for the asymmetry of the ROI. More...
template<typename TPixel >
void	SetImage (TPixel *srcImage, uint srcpitch)
	Set the image on which to perform the tracking. More...
void	SetImage16Bit (ushort *srcImage, uint srcpitch)
	Set an image with 16 bit type. More...
void	SetImage8Bit (uchar *srcImage, uint srcpitch)
	Set an image with 8 bit type. More...
void	SetImageFloat (float *srcImage)
	Set an image with float type. More...
void	SaveImage (const char *filename)
	Save the tracker's image to a jpg file. More...
vector2f	ComputeMeanAndCOM (float bgcorrection=0.0f)
	Calculate the center of mass of the image. More...
void	ComputeRadialProfile (float *dst, int radialSteps, int angularSteps, float minradius, float maxradius, vector2f center, bool crp, bool *boundaryHit=0, bool normalize=true)
	Wrapper to compute a radial profile. More...
void	ComputeQuadrantProfile (scalar_t *dst, int radialSteps, int angularSteps, int quadrant, float minRadius, float maxRadius, vector2f center, float *radialWeights=0)
	Wrapper to compute a radial profile. More...
float	ComputeZ (vector2f center, int angularSteps, int zlutIndex, bool *boundaryHit=0, float *profile=0, float *cmpprof=0, bool normalizeProfile=true)
	Helper function to calculate the Z position. More...
void	FourierTransform2D ()
	Calculate a 2D fourier transform of the tracker's image. More...
void	FourierRadialProfile (float *dst, int radialSteps, int angularSteps, float minradius, float maxradius)
	Calculate the radial profile based on the fourier transform. More...
void	Normalize (float *image=0)
	Normalize an image. More...
void	SetRadialZLUT (float *data, int planes, int res, int num_zluts, float minradius, float maxradius, bool copyMemory, bool useCorrelation)
	Tell the tracker where in memory the LUT is located. More...
void	SetRadialWeights (float *radweights)
	Set the radial weights to be used for profile comparisons. More...
void	CalculateErrorCurve (double *errorcurve_dest, float *profile, float *zlut_sel)
	Calculate the error curve from a specific profile and LUT. More...
void	CalculateInterpolatedZLUTProfile (float *profile_dest, float z, int zlutIndex)
	Calculates an intermediate ZLUT profile between two planes. More...
float	CalculateErrorFlag (double *prof1, double *prof2)
	WIP. Calculate the error flag between two profiles. More...
float	LUTProfileCompare (float *profile, int zlutIndex, float *cmpProf, LUTProfileMaxComputeMode maxPosMethod, float *fitcurve=0, int *maxPos=0, int frameNum=0)
	Compare a profile to a LUT and calculate the optimum Z value for it. More...
float	LUTProfileCompareAdjustedWeights (float *rprof, int zlutIndex, float z_estim)
	Compare a profile to a LUT and calculate the optimum Z value for it using an initial estimate. More...
float *	GetDebugImage ()
	Get the debug image. More...
void	ApplyOffsetGain (float *offset, float *gain, float offsetFactor, float gainFactor)
	Scale the input image with the background calibration values. More...
void	AllocateQIFFTs (int nsteps)
	Initializes the fft handles. More...
vector3f	QuadrantAlign (vector3f initial, int beadIndex, int angularStepsPerQuadrant, bool &boundaryHit)
	Perform the quadrant align algorithm. More...

Public Attributes
int	width
	ROI width. More...
int	height
	ROI height. More...
int	xcorw
	Cross correlation profile length. More...
int	trackerID
	ID of this tracker (= ID of thread it runs in). More...
float *	srcImage
	Memory region that holds the image on which tracking is to be performed. More...
float *	debugImage
	Memory in which an intermediate image can optionally be place and retrieved. More...
float	mean
	Mean intensity of the ROI. Calculated in ComputeMeanAndCOM. More...
float	stdev
	Standard deviation of values in the ROI. Calculated in ComputeMeanAndCOM. More...
float *	zluts
	Pointer to the first data in the ZLUTs. More...
bool	zlut_memoryOwner
	Flag indicating if this instance is the owner of the zluts memory, or is it external. False in all normal operation. More...
int	zlut_planes
	Number of planes per ZLUT. More...
int	zlut_res
	ZLUT resolution = number of radial steps in a plane. More...
int	zlut_count
	Number of ZLUTs (= # beads) available. More...
float	zlut_minradius
	Minimum radius in pixels to start sampling for a ZLUT profile. More...
float	zlut_maxradius
	Maximum radius in pixels of the ZLUT profile sampling circle. More...
std::vector< float >	zlut_radialweights
	Vector with the radialweights used by the error curve calculation. More...
kissfft< scalar_t > *	qa_fft_forward
	Handle to forward FFT kissfft instance for quadrant align. More...
kissfft< scalar_t > *	qa_fft_backward
	Handle to backward FFT kissfft instance for quadrant align. More...
bool	testRun
	Flag to enable running a test run. More...
XCor1DBuffer *	xcorBuffer
	The handle from which to perform cross correlation calculations. More...
std::vector< vector2f >	quadrantDirs
	Vector with the sampling points for a single quadrant (cos & sin pairs). More...
int	qi_radialsteps
	Number of radialsteps in the QI calculations. More...
kissfft< scalar_t > *	qi_fft_forward
	Handle to forward FFT kissfft instance for QI. More...
kissfft< scalar_t > *	qi_fft_backward
	Handle to backward FFT kissfft instance for QI. More...
FFT2D *	fft2d
	Instance of FFT2D to perform 2D FFTs. More...

Detailed Description

Class with all CPU algorithm implementations.

Has one ROI in its memory (srcImage) and the correct member functions are called from its parent's QueuedCPUTracker::ProcessJob. All functions work on the image in memory. Many of the settings that appear here come from QTrkComputedConfig and are further explained there.

Definition at line 38 of file cpu_tracker.h.

Member Enumeration Documentation

LUTProfileMaxComputeMode

enum CPUTracker::LUTProfileMaxComputeMode

Settings to compute the Z coordinate from the error curve.

Enumerator
LUTProfMaxQuadraticFit	Default. Use a 7-point quadratic fit around the error curve's maximum.
LUTProfMaxSplineFit	Compute a spline fit.
LUTProfMaxSimpleInterp	Unimplemented.

Definition at line 378 of file cpu_tracker.h.

                                  { 
        LUTProfMaxQuadraticFit,         
        LUTProfMaxSplineFit,            
        LUTProfMaxSimpleInterp          
    };

Constructor & Destructor Documentation

CPUTracker()

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

Create an instance of CPUTracker.

Parameters

[in]	w	Image width.
[in]	h	Image height.
[in]	xcorwindow	Cross correlation window size.
[in]	testRun	Flag for test run. See QTrkSettings::testRun.

Definition at line 39 of file cpu_tracker.cpp.

    : width(w), height(h), xcorw(xcorwindow), testRun(testMode)
{
    trackerID = 0;

    xcorBuffer = 0;
    
    mean=0.0f;
    srcImage = new float [w*h];
    debugImage = new float [w*h];
    std::fill(srcImage, srcImage+w*h, 0.0f);
    std::fill(debugImage, debugImage+w*h, 0.0f);

    zluts = 0;
    zlut_planes = zlut_res = zlut_count = 0;
    zlut_minradius = zlut_maxradius = 0.0f;

    qa_fft_forward = qa_fft_backward = 0;

    qi_radialsteps = 0;
    qi_fft_forward = qi_fft_backward = 0;

    fft2d=0;
}

~CPUTracker()

CPUTracker::~CPUTracker

(

)

Destroy a CPUTracker instance.

Definition at line 64 of file cpu_tracker.cpp.

{
    delete[] srcImage;
    delete[] debugImage;
    if (zluts && zlut_memoryOwner) 
        delete[] zluts;

    if (xcorBuffer)
        delete xcorBuffer;

    if (qi_fft_forward) { 
        delete qi_fft_forward;
        delete qi_fft_backward;
    }

    if (qa_fft_backward) {
        delete qa_fft_backward;
        delete qa_fft_forward;
    }

    if (fft2d)
        delete fft2d;
}

Member Function Documentation

AllocateQIFFTs()

void CPUTracker::AllocateQIFFTs

(

int

nsteps

)

Initializes the fft handles.

Readies qi_fft_forward and qi_fft_backward for use.

Parameters

[in]

nsteps

The number of radial steps per profile. FFT length is 2 * nsteps.

Definition at line 247 of file cpu_tracker.cpp.

{
    if(!qi_fft_forward || qi_radialsteps != nr) {
        if(qi_fft_forward) {
            delete qi_fft_forward;
            delete qi_fft_backward;
        }
        qi_radialsteps = nr;
        qi_fft_forward = new kissfft<scalar_t>(nr*2,false);
        qi_fft_backward = new kissfft<scalar_t>(nr*2,true);
    }
}

ApplyOffsetGain()

void CPUTracker::ApplyOffsetGain	(	float *	offset,
		float *	gain,
		float	offsetFactor,
		float	gainFactor
	)

Scale the input image with the background calibration values.

See QueuedTracker::SetPixelCalibrationImages() for more information.

Parameters

[in]	offset	Array with per-pixel offsets.
[in]	gain	Array with per-pixel gains.
[in]	offsetFactor	Factor by which to scale the offsets.
[in]	gainFactor	Factor by which to scale the gains.

Definition at line 95 of file cpu_tracker.cpp.

{
    if (offset && !gain) {
        for (int i=0;i<width*height;i++)
            srcImage[i] = srcImage[i]+offset[i]*offsetFactor;
    }
    if (gain && !offset) {
        for (int i=0;i<width*height;i++)
            srcImage[i] = srcImage[i]*gain[i]*gainFactor;
    }
    if (gain && offset)  {
        for (int i=0;i<width*height;i++)
            srcImage[i] = (srcImage[i]+offset[i]*offsetFactor)*gain[i]*gainFactor;
    }
}

CalculateErrorCurve()

void CPUTracker::CalculateErrorCurve	(	double *	errorcurve_dest,
		float *	profile,
		float *	zlut_sel
	)

Calculate the error curve from a specific profile and LUT.

Parameters

[out]	errorcurve_dest	Pre-allocated output array of size zlut_planes which will be filled with the error curve.
[in]	profile	Array with the profile to compare against the LUT.
[in]	zlut_sel	Pointer to the first element of the selected LUT.

Definition at line 776 of file cpu_tracker.cpp.

{
#if 0
    float* zlut_norm = ALLOCA_ARRAY(float, zlut_res);
    float* prof_norm = ALLOCA_ARRAY(float, zlut_res);
    for (int r=0;r<zlut_res;r++)
        prof_norm[r] = rprof[r] * zlut_radialweights[r];
    NormalizeRadialProfile(prof_norm, zlut_res);

    for (int k=0;k<zlut_planes;k++) {
        double diffsum = 0.0f;

        if (zlut_radialweights.empty()) {
            for (int r=0;r<zlut_res;r++)  zlut_norm[r]=zlut_sel[k*zlut_res+r];
        } else {
            for (int r=0;r<zlut_res;r++) {
                zlut_norm[r] = zlut_sel[k*zlut_res+r] * zlut_radialweights[r];
            }
        }
        NormalizeRadialProfile(zlut_norm, zlut_res);

        for (int r = 0; r<zlut_res;r++) {
            double d = prof_norm[r]-zlut_norm[r];
            d = -d*d;
            diffsum += d;
        }
        rprof_diff[k] = diffsum;
    }
#else
    for (int k=0;k<zlut_planes;k++) {
        double diffsum = 0.0f;
        for (int r = 0; r<zlut_res;r++) {
            double d = profile[r]-zlut_sel[k*zlut_res+r];
            if (!zlut_radialweights.empty())
                d*=zlut_radialweights[r];
            d = -d*d;
            diffsum += d;
        }
        errorcurve_dest[k] = diffsum;
    }
#endif
}

CalculateErrorFlag()

float CPUTracker::CalculateErrorFlag	(	double *	prof1,
		double *	prof2
	)

WIP. Calculate the error flag between two profiles.

Parameters

[in]	prof1	The first profile. Typically the measured profile.
[in]	prof2	The second profile. Typically a profile from the LUT.

Returns

The error value between the two profiles. If 'high', tracking was unreliable.

Definition at line 833 of file cpu_tracker.cpp.

{
    //      sum(abs(expectedCurve-errorCurve))/size(errorCurve,2)
    float flag = 0.0f;
    for(int ii = 0; ii < zlut_res; ii++){
        flag += std::abs(prof1[ii]-prof2[ii]);
    }   
    return flag/zlut_res;
}

CalculateInterpolatedZLUTProfile()

void CPUTracker::CalculateInterpolatedZLUTProfile	(	float *	profile_dest,
		float	z,
		int	zlutIndex
	)

Calculates an intermediate ZLUT profile between two planes.

At each radial step, the value is linearly interpolated based on the two nearest saved profiles. Use to obtain estimated profiles at non-integer z values.

Parameters

[out]	profile_dest	Pre-allocated output array of size zlut_res which will be filled with the error curve.
[in]	z	The height at which to determine the expected profile.
[in]	zlutIndex	The index of the zlut/bead.

Definition at line 819 of file cpu_tracker.cpp.

{
    float* zlut_sel = GetRadialZLUT(zlutIndex);
    for (int r = 0; r<zlut_res;r++) {
        if ( z < 0 )
            profile_dest[r] = 0;
        else if( ((int)z+1) < zlut_planes ) // Index out of bounds if z > LUT resolution
            profile_dest[r] = Lerp(zlut_sel[(int)z*zlut_res+r],zlut_sel[((int)z+1)*zlut_res+r],z-(int)z);
        else
            profile_dest[r] = zlut_sel[(zlut_planes-1)*zlut_res+r];
    }
    NormalizeRadialProfile(profile_dest,zlut_res);
}

CheckBoundaries()

bool CPUTracker::CheckBoundaries	(	vector2f	center,
		float	radius
	)

Test to see if a circle extends outside of image boundaries.

Parameters

[in]	center	The center of the circle to test.
[in]	radius	The radius of the circle.

Returns

False if circle remains inside image boundaries.

Definition at line 157 of file cpu_tracker.cpp.

{
    return center.x + radius >= width ||
        center.x - radius < 0 ||
        center.y + radius >= height ||
        center.y - radius < 0;
}

Compute2DGaussianMLE()

CPUTracker::Gauss2DResult CPUTracker::Compute2DGaussianMLE	(	vector2f	initial,
		int	iterations,
		float	sigma
	)

Calculate a 2D Gaussian fit on the image.

Parameters

[in]	initial	The initial position from which to start the algorithm.
[in]	iterations	The number of iterations to perform.
[in]	sigma	Expected standard deviation of the gaussian fit.

Returns

The fitting results.

Definition at line 497 of file cpu_tracker.cpp.

{
    vector2f pos = initial;
    float I0 = mean*0.5f*width*height;
    float bg = mean*0.5f;

    const float _1oSq2Sigma = 1.0f / (sqrtf(2) * sigma);
    const float _1oSq2PiSigma = 1.0f / (sqrtf(2*3.14159265359f) * sigma);
    const float _1oSq2PiSigma3 = 1.0f / (sqrtf(2*3.14159265359f) * sigma*sigma*sigma);

    for (int i=0;i<iterations;i++)
    {
        double dL_dx = 0.0; 
        double dL_dy = 0.0; 
        double dL_dI0 = 0.0;
        double dL_dIbg = 0.0;
        double dL2_dx = 0.0;
        double dL2_dy = 0.0;
        double dL2_dI0 = 0.0;
        double dL2_dIbg = 0.0;

        double mu_sum = 0.0;
                
        for (int y=0;y<height;y++)
        {
            for (int x=0;x<width;x++)
            {
                float Xexp0 = (x-pos.x + .5f) * _1oSq2Sigma;
                float Yexp0 = (y-pos.y + .5f) * _1oSq2Sigma;
        
                float Xexp1 = (x-pos.x - .5f) * _1oSq2Sigma;
                float Yexp1 = (y-pos.y - .5f) * _1oSq2Sigma;
                
                float DeltaX = 0.5f * erf(Xexp0) - 0.5f * erf(Xexp1);
                float DeltaY = 0.5f * erf(Yexp0) - 0.5f * erf(Yexp1);
                float mu = bg + I0 * DeltaX * DeltaY;
                
                float dmu_dx = I0*_1oSq2PiSigma * ( expf(-Xexp1*Xexp1) - expf(-Xexp0*Xexp0)) * DeltaY;

                float dmu_dy = I0*_1oSq2PiSigma * ( expf(-Yexp1*Yexp1) - expf(-Yexp0*Yexp0)) * DeltaX;
                float dmu_dI0 = DeltaX*DeltaY;
                float dmu_dIbg = 1;
        
                float smp = GetPixel(x,y);
                float f = smp / mu - 1;
                dL_dx += dmu_dx * f;
                dL_dy += dmu_dy * f;
                dL_dI0 += dmu_dI0 * f;
                dL_dIbg += dmu_dIbg * f;

                float d2mu_dx = I0*_1oSq2PiSigma3 * ( (x - pos.x - .5f) * expf (-Xexp1*Xexp1) - (x - pos.x + .5f) * expf(-Xexp0*Xexp0) ) * DeltaY;
                float d2mu_dy = I0*_1oSq2PiSigma3 * ( (y - pos.y - .5f) * expf (-Yexp1*Yexp1) - (y - pos.y + .5f) * expf(-Yexp0*Yexp0) ) * DeltaX;
                dL2_dx += d2mu_dx * f - dmu_dx*dmu_dx * smp / (mu*mu);
                dL2_dy += d2mu_dy * f - dmu_dy*dmu_dy * smp / (mu*mu);
                dL2_dI0 += -dmu_dI0*dmu_dI0 * smp / (mu*mu);
                dL2_dIbg += -smp / (mu*mu);

                mu_sum += mu;
            }
        }

        double mean_mu = mu_sum / (width*height);

        pos.x -= dL_dx / dL2_dx;
        pos.y -= dL_dy / dL2_dy;
        I0 -= dL_dI0 / dL2_dI0;
        bg -= dL_dIbg / dL2_dIbg;
    }

    Gauss2DResult r;
    r.pos = pos;
    r.I0 = I0;
    r.bg = bg;

    return r;
}

ComputeAsymmetry()

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

Find a measure for the asymmetry of the ROI.

The center of mass of each radial line (see further explanation of the QI algorithm) is calculated. The standard deviation of these COMs is returned as a measure of the asymmetry.

See QTrkComputedConfig for more information on the settings.

Parameters

[in]	center	The center of the sampling area.
[in]	radialSteps	The number of radial steps per spoke.
[in]	angularSteps	The number of angular steps to take.
[in]	minRadius	Minimum sampling radius in pixels.
[in]	maxRadius	Maximum sampling radius in pixels.
[in]	dstAngProf	Optional. Pre-allocated array to hold the calculated COMs.

Returns

Standard deviation of the radial spokes' centers of mass.

Definition at line 575 of file cpu_tracker.cpp.

{
    vector2f* radialDirs = (vector2f*)ALLOCA(sizeof(vector2f)*angularSteps);
    for (int j=0;j<angularSteps;j++) {
        float ang = 2*3.141593f*j/(float)angularSteps;
        radialDirs[j] = vector2f (cosf(ang), sinf(ang));
    }

    // profile along single radial direction
    float* rline = (float*)ALLOCA(sizeof(float)*radialSteps);

    if (!dstAngProf) 
        dstAngProf = (float*)ALLOCA(sizeof(float)*radialSteps);

    float rstep = (maxRadius-minRadius) / radialSteps;

    for (int a=0;a<angularSteps;a++) {
        // fill rline
        // angularProfile[a] = rline COM 
        
        float r = minRadius;
        for (int i=0;i<radialSteps;i++) {
            float x = center.x + radialDirs[a].x * r;
            float y = center.y + radialDirs[a].y * r;
            rline[i] = Interpolate(srcImage,width,height, x,y);
            r += rstep;
        }

        // Compute rline COM
        dstAngProf[a] = ComputeBgCorrectedCOM1D(rline, radialSteps, 1.0f);
    }
    
    float stdev = ComputeStdDev(dstAngProf, angularSteps);
    return stdev;
}

ComputeMeanAndCOM()

vector2f CPUTracker::ComputeMeanAndCOM

(

float

bgcorrection = 0.0f

)

Calculate the center of mass of the image.

Also calculates and saves the mean and standard deviation of the image. Always performed as first step in localization for a first guess.

The moments are calculated as absolute values around the mean. An optional background correction is possible by only taking into account outliers.

Parameters

[in]

bgcorrection

Factor of standard deviation not to take into account for the center of mass calculation. Default is 0.

Definition at line 661 of file cpu_tracker.cpp.

{
    double sum=0, sum2=0;
    double momentX=0;
    double momentY=0;
    // Order is important 
    //COM: x:53.959133, y:53.958984
    //COM: x:53.958984, y:53.959133 

    for (int y=0;y<height;y++) {
        for (int x=0;x<width;x++)  {
            float v = GetPixel(x,y);
            sum += v;
            sum2 += v*v;
        }
    }

    float invN = 1.0f/(width*height);
    mean = sum * invN;
    stdev = sqrtf(sum2 * invN - mean * mean);
    sum = 0.0f;

    /*float *ymom = ALLOCA_ARRAY(float, width);
    float *xmom = ALLOCA_ARRAY(float, height);

    for(int x=0;x<width;x++)
        ymom[x]=0;
    for(int y=0;y<height;y++)
        xmom[y]=0;*/

    // Comments: Gaussian mask not currently used
    //vector2f centre = vector2f((float)width/2,(float)height/2);
    //float sigma = (float)width/4; // Suppress to 5% at edges
    //float devi = 1/(2*sigma*sigma);
    
    for(int x=0;x<width;x++) {
        for (int y=0;y<height;y++) {
            float v = GetPixel(x,y);
            
            //float xabs = (x-centre.x); float yabs = (y-centre.y);
            //float gaussfact = 1.0f;//expf(- xabs*xabs*devi - yabs*yabs*devi);
            
            v =  std::max(0.0f, fabs(v-mean)-bgcorrection*stdev);
            sum += v;
    //      xmom[y] += x*v;
    //      ymom[x] += y*v;
            momentX += x*(double)v;
            momentY += y*(double)v;
        }
    }
    
    vector2f com;
    com.x = (float)( momentX / sum );
    com.y = (float)( momentY / sum );
    return com;
}

ComputeQI()

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

Execute the quadrant interpolation algorithm.

See [3] for algorithm details. See QTrkComputedConfig for more information on the settings.

Parameters

[in]	initial	The initial position from which to start the algorithm.
[in]	iterations	The number of iterations to perform.
[in]	radialSteps	Amount of radial steps per quadrant profile.
[in]	angularStepsPerQuadrant	Amount of angular sampling steps per profile.
[in]	angStepIterationFactor	Angular step iteration factor, see QTrkSettings::qi_angstep_factor.
[in]	minRadius	Minimum radius of the sampling area in pixels.
[in]	maxRadius	Maximum radius of the sampling area in pixels.
[in,out]	boundaryHit	Pre-allocated bool that will be set to true if the image boundaries have been exceeded during sampling.
[in]	radialweights	Array of the radial weights to use.

Returns

The resulting position with applied offsets.

Definition at line 260 of file cpu_tracker.cpp.

{
    int nr=radialSteps;
#ifdef _DEBUG
    std::copy(srcImage, srcImage+width*height, debugImage);
    maxImageValue = *std::max_element(srcImage,srcImage+width*height);
#endif

    if (angularStepsPerQ != quadrantDirs.size()) {
        quadrantDirs.resize(angularStepsPerQ);
        for (int j=0;j<angularStepsPerQ;j++) {
            float ang = 0.5*3.141593f*(j+0.5f)/(float)angularStepsPerQ;
            quadrantDirs[j] = vector2f( cosf(ang), sinf(ang) );
        }
    }

    AllocateQIFFTs(radialSteps);

    scalar_t* buf = (scalar_t*)ALLOCA(sizeof(scalar_t)*nr*4);
    scalar_t* q0=buf, *q1=buf+nr, *q2=buf+nr*2, *q3=buf+nr*3;
    complex_t* concat0 = (complex_t*)ALLOCA(sizeof(complex_t)*nr*2);
    complex_t* concat1 = concat0 + nr;

    vector2f center = initial;

    float pixelsPerProfLen = (maxRadius-minRadius)/radialSteps;
    boundaryHit = false;

    float angsteps = angularStepsPerQ / powf(angStepIterationFactor, iterations);
    
    for (int k=0;k<iterations;k++){
        // check bounds
        boundaryHit = CheckBoundaries(center, maxRadius);

        for (int q=0;q<4;q++) {
            ComputeQuadrantProfile(buf+q*nr, nr, angsteps, q, minRadius, maxRadius, center, radialWeights);
            //NormalizeRadialProfile(buf+q*nr, nr);
        }
#ifdef QI_DEBUG
        cmp_cpu_qi_prof.assign (buf,buf+4*nr);
#endif

        // Build Ix = [ qL(-r)  qR(r) ]
        // qL = q1 + q2   (concat0)
        // qR = q0 + q3   (concat1)
        for(int r=0;r<nr;r++) {
            concat0[nr-r-1] = q1[r]+q2[r];
            concat1[r] = q0[r]+q3[r];
        }

        scalar_t offsetX = QI_ComputeOffset(concat0, nr, 0);

        // Build Iy = [ qB(-r)  qT(r) ]
        // qT = q0 + q1
        // qB = q2 + q3
        for(int r=0;r<nr;r++) {
            concat1[r] = q0[r]+q1[r];
            concat0[nr-r-1] = q2[r]+q3[r];
        }
        
        scalar_t offsetY = QI_ComputeOffset(concat0, nr, 1);

#ifdef QI_DBG_EXPORT
        std::copy(concat0, concat0+nr*2,tmp.begin()+nr*2);
        WriteComplexImageAsCSV("cpuprofxy.txt", &tmp[0], nr, 4);
        dbgprintf("[%d] OffsetX: %f, OffsetY: %f\n", k, offsetX, offsetY);
#endif

        center.x += offsetX * pixelsPerProfLen;
        center.y += offsetY * pixelsPerProfLen;

        angsteps *= angStepIterationFactor;
    }

    return center;
}

ComputeQuadrantProfile()

void CPUTracker::ComputeQuadrantProfile	(	scalar_t *	dst,
		int	radialSteps,
		int	angularSteps,
		int	quadrant,
		float	minRadius,
		float	maxRadius,
		vector2f	center,
		float *	radialWeights = 0
	)

Wrapper to compute a radial profile.

Parameters

[in,out]	dst	Pre-allocated float array of size radialSteps to return the profile.
[in]	radialSteps	Number of radial steps in the profile.
[in]	angularSteps	Number of angular steps in the profile.
[in]	quadrant	Quadrant number. 0-3, 1 is upper right, progresses in counterclockwise order.
[in]	minRadius	Minimum sampling circle radius in pixels.
[in]	maxRadius	Maximum sampling circle radius in pixels.
[in]	center	The center around which to sample.
[in]	radialWeights	Optional. Array of radial weights by which to weigh the profile.

Definition at line 613 of file cpu_tracker.cpp.

{
    const int qmat[] = {
        1, 1,
        -1, 1,
        -1, -1,
        1, -1 };
    int mx = qmat[2*quadrant+0];
    int my = qmat[2*quadrant+1];

    if (angularSteps < MIN_RADPROFILE_SMP_COUNT)
        angularSteps = MIN_RADPROFILE_SMP_COUNT;

    for (int i=0;i<radialSteps;i++)
        dst[i]=0.0f;

    scalar_t total = 0.0f;
    scalar_t rstep = (maxRadius - minRadius) / radialSteps;
    for (int i=0;i<radialSteps; i++) {
        scalar_t sum = 0.0f;
        scalar_t r = minRadius + rstep * i;

        int nPixels = 0;
        scalar_t angstepf = (scalar_t) quadrantDirs.size() / angularSteps;
        for (int a=0;a<angularSteps;a++) {
            int i = (int)angstepf * a;
            scalar_t x = center.x + mx*quadrantDirs[i].x * r;
            scalar_t y = center.y + my*quadrantDirs[i].y * r;

            bool outside;
            scalar_t v = Interpolate(srcImage,width,height, x,y, &outside);
            if (!outside) {
                sum += v;
                nPixels++;
                MARKPIXELI(x,y);
            }
        }

        dst[i] = nPixels>=MIN_RADPROFILE_SMP_COUNT ? sum/nPixels : mean;
        if (radialWeights) dst[i] *= radialWeights[i];
        total += dst[i];
    }
    
    //WriteImageAsCSV(SPrintf("D:\\TestImages\\qprof%d.txt", quadrant).c_str(), dst, 1, radialSteps);   
    //WriteArrayAsCSVRow("D:\\TestImages\\radialWeights.csv",radialWeights,radialSteps,false);
}

ComputeRadialProfile()

void CPUTracker::ComputeRadialProfile	(	float *	dst,
		int	radialSteps,
		int	angularSteps,
		float	minradius,
		float	maxradius,
		vector2f	center,
		bool	crp,
		bool *	boundaryHit = 0,
		bool	normalize = true
	)

Wrapper to compute a radial profile.

Parameters

[in,out]	dst	Pre-allocated float array of size radialSteps to return the profile.
[in]	radialSteps	Number of radial steps in the profile.
[in]	angularSteps	Number of angular steps in the profile.
[in]	minradius	Minimum sampling circle radius in pixels.
[in]	maxradius	Maximum sampling circle radius in pixels.
[in]	center	The center around which to sample.
[in]	crp	No clue, good luck Josko!
[in,out]	boundaryHit	Optional. Bool set to indicate whether the image boundary has been hit during sampling.
[in]	normalize	Optional. Normalize the radial profile?

Definition at line 727 of file cpu_tracker.cpp.

{
    bool boundaryHit = CheckBoundaries(center, maxradius);
    if (pBoundaryHit) *pBoundaryHit = boundaryHit;

    ImageData imgData (srcImage, width,height);
    if (crp) {
        ComputeCRP(dst, radialSteps, angularSteps, minradius, maxradius, center, &imgData, mean);
    }
    else {
        ::ComputeRadialProfile(dst, radialSteps, angularSteps, minradius, maxradius, center, &imgData, mean, normalize);
    }
}

ComputeXCorInterpolated()

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

Compute the cross correlation offsets and resulting position.

See https://en.wikipedia.org/wiki/Cross-correlation for algorithm details.

Parameters

[in]	initial	The initial position from which to start the algorithm.
[in]	iterations	The number of iterations to perform.
[in]	profileWidth	Width of the profiles to correlate.
[in,out]	boundaryHit	Pre-allocated bool that will be set to true if the image boundaries have been exceeded during sampling.

Returns

The resulting position with applied offsets.

Definition at line 165 of file cpu_tracker.cpp.

{
    // extract the image
    vector2f pos = initial;

    if (!xcorBuffer)
        xcorBuffer = new XCor1DBuffer(xcorw);

    if (xcorw < profileWidth)
        profileWidth = xcorw;

#ifdef _DEBUG
    std::copy(srcImage, srcImage+width*height, debugImage);
    maxImageValue = *std::max_element(srcImage,srcImage+width*height);
#endif

    if (xcorw > width || xcorw > height) {
        boundaryHit = true;
        return initial;
    }

    complex_t* prof = (complex_t*)ALLOCA(sizeof(complex_t)*xcorw);
    complex_t* prof_rev = (complex_t*)ALLOCA(sizeof(complex_t)*xcorw);
    scalar_t* prof_autocor = (scalar_t*)ALLOCA(sizeof(scalar_t)*xcorw);

    boundaryHit = false;
    for (int k=0;k<iterations;k++) {
        boundaryHit = CheckBoundaries(pos, XCorScale*xcorw/2);

        float xmin = pos.x - XCorScale * xcorw/2;
        float ymin = pos.y - XCorScale * xcorw/2;

        // generate X position xcor array (summing over y range)
        for (int x=0;x<xcorw;x++) {
            scalar_t s = 0.0f;
            int n=0;
            for (int y=0;y<profileWidth;y++) {
                scalar_t xp = x * XCorScale + xmin;
                scalar_t yp = pos.y + XCorScale * (y - profileWidth/2);
                bool outside;
                s += Interpolate(srcImage, width, height, xp, yp, &outside);
                n += outside?0:1;
                MARKPIXELI(xp, yp);
            }
            if (n >0) s/=n;
            else s=0;
            prof [x] = s;
            prof_rev [xcorw-x-1] = s;
        }

        xcorBuffer->XCorFFTHelper(prof, prof_rev, prof_autocor);
        scalar_t offsetX = ComputeMaxInterp<scalar_t,QI_LSQFIT_NWEIGHTS>::Compute(prof_autocor, xcorw, QIWeights) - (scalar_t)xcorw/2;

        // generate Y position xcor array (summing over x range)
        for (int y=0;y<xcorw;y++) {
            scalar_t s = 0.0f; 
            int n=0;
            for (int x=0;x<profileWidth;x++) {
                scalar_t xp = pos.x + XCorScale * (x - profileWidth/2);
                scalar_t yp = y * XCorScale + ymin;
                bool outside;
                s += Interpolate(srcImage,width,height, xp, yp, &outside);
                n += outside?0:1;
                MARKPIXELI(xp,yp);
            }
            if (n >0) s/=n;
            else s=0;
            prof[y] = s;
            prof_rev [xcorw-y-1] =s;
        }

        xcorBuffer->XCorFFTHelper(prof,prof_rev, prof_autocor);
        //WriteImageAsCSV("xcorautoconv.txt",&xcorBuffer->Y_result[0],xcorBuffer->Y_result.size(),1);
        scalar_t offsetY = ComputeMaxInterp<scalar_t, QI_LSQFIT_NWEIGHTS>::Compute(prof_autocor, xcorw, QIWeights) - (scalar_t)xcorw/2;

        pos.x += (offsetX - 1) * XCorScale * 0.5f;
        pos.y += (offsetY - 1) * XCorScale * 0.5f;
    }

    return pos;
}

ComputeZ()

float CPUTracker::ComputeZ ( vector2f  center, int  angularSteps, int  zlutIndex, bool *  boundaryHit = 0, float *  profile = 0, float *  cmpprof = 0, bool  normalizeProfile = true  )

inline

Helper function to calculate the Z position.

Calculates the radial profile with CPUTracker::ComputeRadialProfile and calls CPUTracker::LUTProfileCompare to compare it to the LUT.

Parameters

[in]	center	The center around which to sample.
[in]	angularSteps	Number of angular steps in the profile.
[in]	zlutIndex	Index of the ZLUT/bead number.
[in,out]	boundaryHit	Optional. Bool set to indicate whether the image boundary has been hit during sampling.
[in,out]	profile	Optional. Pre-allocated array to retrieve the profile if so desired.
[in,out]	cmpprof	Optional. Pre-allocated array to retrieve the error curve if so desired..
[in]	normalizeProfile	Optional. Normalize the profile? Default is true.

Definition at line 319 of file cpu_tracker.h.

    {
        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);
    }

FourierRadialProfile()

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

Calculate the radial profile based on the fourier transform.

Setting is available, not ever really used.

Parameters

[in,out]	dst	Pre-allocated float array of size radialSteps to return the profile.
[in]	radialSteps	Number of radial steps in the profile.
[in]	angularSteps	Number of angular steps in the profile.
[in]	minradius	Minimum sampling circle radius in pixels.
[in]	maxradius	Maximum sampling circle radius in pixels.

Definition at line 1060 of file cpu_tracker.cpp.

{
    FourierTransform2D();
    ImageData img(srcImage,width,height);
    ::ComputeRadialProfile(dst, radialSteps, angularSteps, 1, width*0.3f, vector2f(width/2,height/2), &img, 0, false);
    for (int i=0;i<radialSteps;i++) {
        dst[i]=logf(dst[i]);
    }
    //NormalizeRadialProfile(dst, radialSteps);
}

FourierTransform2D()

void CPUTracker::FourierTransform2D

(

)

Calculate a 2D fourier transform of the tracker's image.

See FFT2D for more information.

Definition at line 1016 of file cpu_tracker.cpp.

{
    //kissfft<float> yfft(h, inverse);
    if (!fft2d){
        fft2d = new FFT2D(width, height);
    }

    fft2d->Apply(srcImage);
}

GetDebugImage()

float* CPUTracker::GetDebugImage ( )

inline

Get the debug image.

The debug image can be used to store intermediate results of image processing and retrieved using this function.

Returns

A pointer to the array with the debug image in it.

Definition at line 440 of file cpu_tracker.h.

{ return debugImage; }

GetHeight()

int CPUTracker::GetHeight ( )

inline

Get the height of the image.

Definition at line 124 of file cpu_tracker.h.

GetPixel()

float& CPUTracker::GetPixel ( int  x, int  y  )

inline

Get the image pixel greyscale value at point (x,y).

Definition at line 122 of file cpu_tracker.h.

GetRadialZLUT()

float* CPUTracker::GetRadialZLUT ( int  index )

inline

Get the start of the ZLUT of a specific bead.

Parameters

[in]

index

The bead number for which to get the LUT.

Returns

A pointer to the start of the LUT of this bead.

Definition at line 90 of file cpu_tracker.h.

{ return &zluts[zlut_res*zlut_planes*index]; }

GetWidth()

int CPUTracker::GetWidth ( )

inline

Get the width of the image.

Definition at line 123 of file cpu_tracker.h.

LUTProfileCompare()

float CPUTracker::LUTProfileCompare	(	float *	profile,
		int	zlutIndex,
		float *	cmpProf,
		LUTProfileMaxComputeMode	maxPosMethod,
		float *	fitcurve = 0,
		int *	maxPos = 0,
		int	frameNum = 0
	)

Compare a profile to a LUT and calculate the optimum Z value for it.

Parameters

[in]	profile	The profile to compare.
[in]	zlutIndex	The number of the ZLUT to compare with.
[in]	cmpProf	Pre-allocated array in which to return the error curve. Not used if null pointer is passed.
[in]	maxPosMethod	Which method to use to determine the Z position from the error curve. See LUTProfileMaxComputeMode.
[in]	fitcurve	Optional. Pre-allocated array in which to return the error curve fit. Not used if null pointer is passed.
[in]	maxPos	Optional. Pre-allocated integer in which the position of the error curve's maximum will be put.
[in]	frameNum	Optional. The number of the frame this track job belongs to. Only used in testRun.

Definition at line 843 of file cpu_tracker.cpp.

{
    /* From this function save:
        - Frame number?
        
        - Original image
        - Profile
        - LUT
        - Error curve
        - Final value
    */

    bool freeMem = false;
    std::string outputName;
    float zOffset;

    if(testRun){
        outputName = SPrintf("%s%05d-%05d",GetCurrentOutputPath().c_str(),zlutIndex,frameNum);
        
        if (FILE *file = fopen(SPrintf("%s.jpg",outputName.c_str()).c_str(), "r")) {
            fclose(file);
            printf("File exists\n");
        }
        SaveImage(SPrintf("%s.jpg",outputName.c_str()).c_str());
        //FloatToJPEGFile(SPrintf("%s-prof.jpg",outputName.c_str()).c_str(), rprof, zlut_res, 1);
        WriteArrayAsCSVRow(SPrintf("%s-prof.txt",outputName.c_str()).c_str(), rprof, zlut_res, 0);

        if(!fitcurve){
            fitcurve = new float[zlut_planes];
            freeMem = true;
        }

        /* 
            For testing, we want to introduce a fixed height offset of xx nm, regardless of amount of zlut planes
            Z is in zlut planes here, so calculate back

            int numImgInStack = 1218;
            int numPositions = 1001; // ~10nm/frame
            float range = 10.0f; // total range 25.0 um -> 35.0 um
            float umPerImg = range/numImgInStack; //
            
            int zplanes = 50;
            int startFrame = 400;
        */
        float range = 10.0f; // total range 25.0 um -> 35.0 um
        int numImgInStack = 1218;
        float umPerImg = range/numImgInStack;
        int startFrame = 400;
        int usedFrames = numImgInStack-startFrame;
        float rangeInLUT = usedFrames*umPerImg; // Total um of LUT range
        float umPerPlane = rangeInLUT/zlut_planes;

        float error = 0.000f; // Forced error in um
        zOffset = error/umPerPlane; // um * planes / um
    }

    if (!zluts)
        return 0.0f;

    double* rprof_diff = ALLOCA_ARRAY(double, zlut_planes);
    //WriteImageAsCSV("zlutradprof-cpu.txt", rprof, zlut_res, 1);

    // Now compare the radial profile to the profiles stored in Z
    float* zlut_sel = GetRadialZLUT(zlutIndex);

    if(testRun){
        FloatToJPEGFile(SPrintf("%s-zlut.jpg",outputName.c_str()).c_str(),zlut_sel,zlut_res,zlut_planes);
    }

    CalculateErrorCurve(rprof_diff, rprof, zlut_sel);

    if (cmpProf) {
        //cmpProf->resize(zlut_planes);
        std::copy(rprof_diff, rprof_diff+zlut_planes, cmpProf);
    }

    if (maxPosComputeMode == LUTProfMaxQuadraticFit) {
        if (maxPos) {
            *maxPos = std::max_element(rprof_diff, rprof_diff+zlut_planes) - rprof_diff;
        }

        if (fitcurve) {
            LsqSqQuadFit<double> fit;
            float z = ComputeMaxInterp<double, ZLUT_LSQFIT_NWEIGHTS>::Compute(rprof_diff, zlut_planes, ZLUTWeights_d, &fit);

            int iMax = std::max_element(rprof_diff, rprof_diff+zlut_planes) - rprof_diff;
            for (int i=0;i<zlut_planes;i++)
                fitcurve[i] = fit.compute(i-iMax);
            
            if(testRun) {
                z += zOffset; // For testing purposes

                // Calculate errorcurve of found ZLUT plane
                float* int_prof = ALLOCA_ARRAY(float, zlut_res);
                CalculateInterpolatedZLUTProfile(int_prof, z, zlutIndex);

                double* plane_auto_diff = ALLOCA_ARRAY(double, zlut_planes);
                CalculateErrorCurve(plane_auto_diff, int_prof, zlut_sel);
                /*for (int k=0;k<zlut_planes;k++) {
                    double diffsum = 0.0f;
                    for (int r = 0; r<zlut_res;r++) {
                        double d = int_prof[r] - zlut_sel[k*zlut_res+r];
                        if (!zlut_radialweights.empty())
                            d*=zlut_radialweights[r];
                        d = -d*d;
                        diffsum += d;
                    }
                    plane_auto_diff[k] = diffsum;
                }*/

                //float errorFlag = CalculateErrorFlag(rprof_diff, plane_auto_diff);

                WriteArrayAsCSVRow(SPrintf("%s-int-prof.txt",outputName.c_str()).c_str(), int_prof, zlut_res, 0);

                FILE* f = fopen(SPrintf("%s-out.txt",outputName.c_str()).c_str(),"w+");
                fprintf_s(f,"z: %f\n",z);
                for (int i=0;i<zlut_planes;i++)
                    fprintf_s(f,"%f\t%f\t%f\t%f\t%f\n",rprof_diff[i],fitcurve[i],plane_auto_diff[i],fit.a,fit.vertexForm());
                fclose(f);

                if(freeMem)
                    delete[] fitcurve;
            }
            
            return z;
        } else {
            return ComputeMaxInterp<double, ZLUT_LSQFIT_NWEIGHTS>::Compute(rprof_diff, zlut_planes, ZLUTWeights_d);
        }
    }
    else if (maxPosComputeMode == LUTProfMaxSimpleInterp) {
        assert(0);
        return 0;
    } else 
        return ComputeSplineFitMaxPos(rprof_diff, zlut_planes);
}

LUTProfileCompareAdjustedWeights()

float CPUTracker::LUTProfileCompareAdjustedWeights	(	float *	rprof,
		int	zlutIndex,
		float	z_estim
	)

Compare a profile to a LUT and calculate the optimum Z value for it using an initial estimate.

Warning

Specific use unknown, only called when LT_LocalizeZWeighted is enabled. Probably was used to test the effect of different weights and fit methods.

Parameters

[in]	rprof	The profile to compare.
[in]	zlutIndex	The number of the ZLUT to compare against.
[in]	z_estim	An initial guess for the z position.

Definition at line 979 of file cpu_tracker.cpp.

{
    if (!zluts)
        return 0.0f;
    
    double* rprof_diff = ALLOCA_ARRAY(double, zlut_planes);

    // Now compare the radial profile to the profiles stored in Z
    float* zlut_sel = GetRadialZLUT(zlutIndex);

    int zi = (int)z_estim;
    if (zi > zlut_planes-2) zi = zlut_planes-2;
    float* z0 = &zlut_sel[zi*zlut_res];
    float* z1 = &zlut_sel[(zi+1)*zlut_res];
    double* zlut_weights = ALLOCA_ARRAY(double, zlut_res);
    for (int i=0;i<zlut_res;i++)
        zlut_weights[i] = fabs((double)z1[i]-(double)z0[i]) * ( zlut_minradius + (zlut_maxradius - zlut_minradius) * i/(float)zlut_res );

    for (int k=0;k<zlut_planes;k++) {
        double diffsum = 0.0f;
        for (int r = 0; r<zlut_res;r++) {
            double d = (double)rprof[r]-(double)zlut_sel[k*zlut_res+r];
            d = -d*d;
            d *= zlut_weights[r];
            diffsum += d;
        }
        rprof_diff[k] = diffsum;
    }
    //return ComputeSplineFitMaxPos(rprof_diff,zlut_planes);
    return ComputeMaxInterp<double, ZLUT_LSQFIT_NWEIGHTS>::Compute(rprof_diff, zlut_planes, ZLUTWeights_d);
}

Normalize()

void CPUTracker::Normalize

(

float *

image = 0

)

Normalize an image.

See normalize.

Parameters

[in]

image

Optional. The image to normalize. If a null pointer (default) is given, uses the source image.

Definition at line 719 of file cpu_tracker.cpp.

{
    if (!d) d=srcImage;
    normalize(d, width, height);
}

QI_ComputeOffset()

scalar_t CPUTracker::QI_ComputeOffset	(	complex_t *	qi_profile,
		int	nr,
		int	axisForDebug
	)

QI helper function to calculate the offset from concatenated profiles.

Parameters

[in]	qi_profile	The concatenated profiles' FFT from which to obtain the offset.
[in]	nr	The number of radial steps per quadrant profile.
[in]	axisForDebug	Axis number used to save intermediate data.

Returns

The calculated offset.

Definition at line 351 of file cpu_tracker.cpp.

{
    complex_t* reverse = ALLOCA_ARRAY(complex_t, nr*2);
    complex_t* fft_out = ALLOCA_ARRAY(complex_t, nr*2);
    complex_t* fft_out2 = ALLOCA_ARRAY(complex_t, nr*2);

    for(int x=0;x<nr*2;x++)
        reverse[x] = profile[nr*2-1-x];

    qi_fft_forward->transform(profile, fft_out);
    qi_fft_forward->transform(reverse, fft_out2); // fft_out2 contains fourier-domain version of reverse profile

    // multiply with conjugate
    for(int x=0;x<nr*2;x++)
        fft_out[x] *= conjugate(fft_out2[x]);

    qi_fft_backward->transform(fft_out, fft_out2);

#ifdef QI_DEBUG
    cmp_cpu_qi_fft_out.assign(fft_out2, fft_out2+nr*2);
#endif

    // fft_out2 now contains the autoconvolution
    // convert it to float
    scalar_t* autoconv = ALLOCA_ARRAY(scalar_t, nr*2);
    for(int x=0;x<nr*2;x++)  {
        autoconv[x] = fft_out2[(x+nr)%(nr*2)].real();
    }

    float* profileReal = ALLOCA_ARRAY(float, nr*2);
    for(int x=0;x<nr*2;x++)  {
        profileReal[x] = profile[x].real();
    }
    //WriteArrayAsCSVRow(SPrintf("D:\\TestImages\\AutoconvProf-%d.txt",axisForDebug).c_str()    ,profileReal,nr*2,false);
    //WriteArrayAsCSVRow(SPrintf("D:\\TestImages\\Autoconv-%d.txt",axisForDebug).c_str()        ,autoconv,nr*2,false); // */
    
    scalar_t maxPos = ComputeMaxInterp<scalar_t, QI_LSQFIT_NWEIGHTS>::Compute(autoconv, nr*2, QIWeights);
    return (maxPos - nr) * (3.14159265359f / 4);
}

QuadrantAlign()

vector3f CPUTracker::QuadrantAlign	(	vector3f	initial,
		int	beadIndex,
		int	angularStepsPerQuadrant,
		bool &	boundaryHit
	)

Perform the quadrant align algorithm.

See LT_ZLUTAlign.

Parameters

[in]	initial	The initial, 3 dimensional position.
[in]	beadIndex	The index of the bead/zlut.
[in]	angularStepsPerQuadrant	Number of angular steps per quadrant.
[in]	boundaryHit	Pre-allocated bool to return whether the image boundary has been hit.

Returns

The updated 3D position.

Definition at line 435 of file cpu_tracker.cpp.

{
    float* zlut = GetRadialZLUT(beadIndex);
    int res=zlut_res;
    complex_t* profile = ALLOCA_ARRAY(complex_t, res*2);
    complex_t* concat0 = ALLOCA_ARRAY(complex_t, res*2);
    complex_t* concat1 = concat0 + res;

    memset(concat0, 0, sizeof(complex_t)*res*2);

    int zp0 = clamp( (int)pos.z, 0, zlut_planes - 1);
    int zp1 = clamp( (int)pos.z + 1, 0, zlut_planes - 1);

    float* zlut0 = &zlut[ res * zp0 ];
    float* zlut1 = &zlut[ res * zp1 ];
    float frac = pos.z - (int)pos.z;
    for (int r=0;r<zlut_res;r++) {
    // Interpolate plane
        double zlutValue = Lerp(zlut0[r], zlut1[r], frac);
        concat0[res-r-1] = concat1[r] = zlutValue;
    }
//  WriteComplexImageAsCSV("qa_zlutprof.txt", concat0, res,2);
    qa_fft_forward->transform(concat0, profile);

    scalar_t* buf = ALLOCA_ARRAY(scalar_t, res*4);
    scalar_t* q0=buf, *q1=buf+res, *q2=buf+res*2, *q3=buf+res*3;

    boundaryHit = CheckBoundaries(vector2f(pos.x,pos.y), zlut_maxradius);
    for (int q=0;q<4;q++) {
        scalar_t *quadrantProfile = buf+q*res;
        ComputeQuadrantProfile(quadrantProfile, res, angularStepsPerQuadrant, q, zlut_minradius, zlut_maxradius, vector2f(pos.x,pos.y));
        NormalizeRadialProfile(quadrantProfile, res);
    }
    
    //WriteImageAsCSV("qa_qdr.txt" , buf, res,4);
    
    float pixelsPerProfLen = (zlut_maxradius-zlut_minradius)/zlut_res;
    boundaryHit = false;

    // Build Ix = [ qL(-r)  qR(r) ]
    // qL = q1 + q2   (concat0)
    // qR = q0 + q3   (concat1)
    for(int r=0;r<res;r++) {
        concat0[res-r-1] = q1[r]+q2[r];
        concat1[r] = q0[r]+q3[r];
    }

    scalar_t offsetX = QuadrantAlign_ComputeOffset(concat0, profile, res, 0);

    // Build Iy = [ qB(-r)  qT(r) ]
    // qT = q0 + q1
    // qB = q2 + q3
    for(int r=0;r<res;r++) {
        concat1[r] = q0[r]+q1[r];
        concat0[res-r-1] = q2[r]+q3[r];
    }
        
    scalar_t offsetY = QuadrantAlign_ComputeOffset(concat0, profile, res, 1);

    return vector3f(pos.x + offsetX * pixelsPerProfLen, pos.y + offsetY * pixelsPerProfLen, pos.z);
}

QuadrantAlign_ComputeOffset()

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

QA helper function to calculate the offset from concatenated profiles.

Parameters

[in]	profile	The concatenated profiles' FFT from which to obtain the offset.
[in]	zlut_prof_fft	The ZLUT profiles' FFT from which to obtain the offset.
[in]	nr	The number of radial steps per quadrant profile.
[in]	axisForDebug	Axis number used to save intermediate data.

Returns

The calculated offset.

Definition at line 393 of file cpu_tracker.cpp.

{
    complex_t* fft_out = ALLOCA_ARRAY(complex_t, nr*2);

//  WriteComplexImageAsCSV("qa_profile.txt", profile, nr*2, 1);

    qa_fft_forward->transform(profile, fft_out);

    // multiply with conjugate
    for(int x=0;x<nr*2;x++)
        fft_out[x] *= conjugate(zlut_prof_fft[x]);

    qa_fft_backward->transform(fft_out, profile);

#ifdef QI_DEBUG
    cmp_cpu_qi_fft_out.assign(fft_out2, fft_out2+nr*2);
#endif

    // profile now contains the autoconvolution
    // convert it to float
    scalar_t* autoconv = ALLOCA_ARRAY(scalar_t, nr*2);
    for(int x=0;x<nr*2;x++)  {
        autoconv[x] = profile[(x+nr)%(nr*2)].real();
    }

//  WriteImageAsCSV("qa_autoconv.txt", autoconv, nr*2, 1);

    scalar_t maxPos = ComputeMaxInterp<scalar_t, QI_LSQFIT_NWEIGHTS>::Compute(autoconv, nr*2, QIWeights);
    return (maxPos - nr) * (3.14159265359f / 4);
    
}

SaveImage()

void CPUTracker::SaveImage

(

const char *

filename

)

Save the tracker's image to a jpg file.

Parameters

[in]

filename

Name of the file. Can also be full path.

Definition at line 1011 of file cpu_tracker.cpp.

{
    FloatToJPEGFile(filename, srcImage, width, height);
}

SetImage()

template<typename TPixel >

void CPUTracker::SetImage	(	TPixel *	srcImage,
		uint	srcpitch
	)

Set the image on which to perform the tracking.

Parameters

[in]	srcImage	Array with the image data.
[in]	srcpitch	Width of one row of data in bytes (typically size(dataType)*imageWidth).

Definition at line 476 of file cpu_tracker.h.

{
    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;
}

SetImage16Bit()

void CPUTracker::SetImage16Bit ( ushort *  srcImage, uint  srcpitch  )

inline

Set an image with 16 bit type.

Parameters

[in]	srcImage	Array with the image data.
[in]	srcpitch	Width of one row of data in bytes (typically size(dataType)*imageWidth).

Definition at line 247 of file cpu_tracker.h.

{ SetImage(srcImage, srcpitch); }

SetImage8Bit()

void CPUTracker::SetImage8Bit ( uchar *  srcImage, uint  srcpitch  )

inline

Set an image with 8 bit type.

Parameters

[in]	srcImage	Array with the image data.
[in]	srcpitch	Width of one row of data in bytes (typically size(dataType)*imageWidth).

Definition at line 254 of file cpu_tracker.h.

{ SetImage(srcImage, srcpitch); }

SetImageFloat()

void CPUTracker::SetImageFloat

(

float *

srcImage

)

Set an image with float type.

Parameters

[in]

srcImage

Array with the image data.

Definition at line 88 of file cpu_tracker.cpp.

{
    for (int k=0;k<width*height;k++)
        srcImage[k]=src[k];
    mean=0.0f;
}

SetRadialWeights()

void CPUTracker::SetRadialWeights

(

float *

radweights

)

Set the radial weights to be used for profile comparisons.

Parameters

[in]

radweights

Array with the radial weights.

Definition at line 767 of file cpu_tracker.cpp.

{
    if (radweights) {
        zlut_radialweights.resize(zlut_res);
        std::copy(radweights, radweights+zlut_res, zlut_radialweights.begin());
    } else
        zlut_radialweights.clear();
}

SetRadialZLUT()

void CPUTracker::SetRadialZLUT	(	float *	data,
		int	planes,
		int	res,
		int	num_zluts,
		float	minradius,
		float	maxradius,
		bool	copyMemory,
		bool	useCorrelation
	)

Tell the tracker where in memory the LUT is located.

The LUT is a large contiguous memory area in which all LUTs are saved. See QueuedCPUTracker::zluts.

It is possible to make one CPUTracker instance the owner of the ZLUT memory if no other manager (like QueuedCPUTracker) is used. For this, use the copyMemory parameter.

Parameters

[in]	data	Pointer to the first element of the ZLUT memory.
[in]	planes	The number of planes per lookup table.
[in]	res	Number of radial steps per plane.
[in]	num_zluts	The number of LUTs in memory.
[in]	minradius	Starting sampling distance of the profiles in pixels.
[in]	maxradius	Ending sampling distance of the profiles in pixels.
[in]	copyMemory	Bool to indicate whether the data should be copied to locally managed memory.
[in]	useCorrelation	Deprecated. No effect, only there not to break calls.

Definition at line 741 of file cpu_tracker.cpp.

{
    if (zluts && zlut_memoryOwner)
        delete[] zluts;
    
    if (copyMemory) {
        zluts = new float[planes*res*numLUTs];
        std::copy(data, data+(planes*res*numLUTs), zluts);
    } else
        zluts = data;
    zlut_memoryOwner = copyMemory;
    zlut_planes = planes;
    zlut_res = res;
    zlut_count = numLUTs;
    zlut_minradius = minradius;
    zlut_maxradius = maxradius;

    if (qa_fft_backward) {
        delete qa_fft_backward;
        delete qa_fft_forward;
    }

    qa_fft_forward = new kissfft<scalar_t>(res*2,false);
    qa_fft_backward = new kissfft<scalar_t>(res*2,true);
}

Member Data Documentation

debugImage

float* CPUTracker::debugImage

Memory in which an intermediate image can optionally be place and retrieved.

Definition at line 47 of file cpu_tracker.h.

fft2d

FFT2D* CPUTracker::fft2d

Instance of FFT2D to perform 2D FFTs.

Definition at line 120 of file cpu_tracker.h.

height

int CPUTracker::height

ROI height.

Definition at line 42 of file cpu_tracker.h.

mean

float CPUTracker::mean

Mean intensity of the ROI. Calculated in ComputeMeanAndCOM.

Definition at line 48 of file cpu_tracker.h.

qa_fft_backward

kissfft<scalar_t>* CPUTracker::qa_fft_backward

Handle to backward FFT kissfft instance for quadrant align.

Definition at line 70 of file cpu_tracker.h.

qa_fft_forward

kissfft<scalar_t>* CPUTracker::qa_fft_forward

Handle to forward FFT kissfft instance for quadrant align.

Definition at line 69 of file cpu_tracker.h.

qi_fft_backward

kissfft<scalar_t>* CPUTracker::qi_fft_backward

Handle to backward FFT kissfft instance for QI.

Definition at line 96 of file cpu_tracker.h.

qi_fft_forward

kissfft<scalar_t>* CPUTracker::qi_fft_forward

Handle to forward FFT kissfft instance for QI.

Definition at line 95 of file cpu_tracker.h.

qi_radialsteps

int CPUTracker::qi_radialsteps

Number of radialsteps in the QI calculations.

Definition at line 94 of file cpu_tracker.h.

quadrantDirs

std::vector<vector2f> CPUTracker::quadrantDirs

Vector with the sampling points for a single quadrant (cos & sin pairs).

Definition at line 93 of file cpu_tracker.h.

srcImage

float* CPUTracker::srcImage

Memory region that holds the image on which tracking is to be performed.

Definition at line 46 of file cpu_tracker.h.

stdev

float CPUTracker::stdev

Standard deviation of values in the ROI. Calculated in ComputeMeanAndCOM.

Definition at line 49 of file cpu_tracker.h.

testRun

bool CPUTracker::testRun

Flag to enable running a test run.

A test run outputs a lot of intermediate data into an output path using GetCurrentOutputPath. This can currently be used to analyse the particulars of the ZLUT algorithms.

Output currently happens in LUTProfileCompare. A Matlab GUI has been created to read and analyse the created files.

Warning

No link to LabVIEW, test executables only.

Definition at line 82 of file cpu_tracker.h.

trackerID

int CPUTracker::trackerID

ID of this tracker (= ID of thread it runs in).

Definition at line 44 of file cpu_tracker.h.

width

int CPUTracker::width

ROI width.

Definition at line 41 of file cpu_tracker.h.

xcorBuffer

XCor1DBuffer* CPUTracker::xcorBuffer

The handle from which to perform cross correlation calculations.

Definition at line 92 of file cpu_tracker.h.

xcorw

int CPUTracker::xcorw

Cross correlation profile length.

Definition at line 43 of file cpu_tracker.h.

zlut_count

int CPUTracker::zlut_count

Number of ZLUTs (= # beads) available.

Definition at line 65 of file cpu_tracker.h.

zlut_maxradius

float CPUTracker::zlut_maxradius

Maximum radius in pixels of the ZLUT profile sampling circle.

Definition at line 67 of file cpu_tracker.h.

zlut_memoryOwner

bool CPUTracker::zlut_memoryOwner

Flag indicating if this instance is the owner of the zluts memory, or is it external. False in all normal operation.

Definition at line 62 of file cpu_tracker.h.

zlut_minradius

float CPUTracker::zlut_minradius

Minimum radius in pixels to start sampling for a ZLUT profile.

Definition at line 66 of file cpu_tracker.h.

zlut_planes

int CPUTracker::zlut_planes

Number of planes per ZLUT.

Definition at line 63 of file cpu_tracker.h.

zlut_radialweights

std::vector<float> CPUTracker::zlut_radialweights

Vector with the radialweights used by the error curve calculation.

Definition at line 68 of file cpu_tracker.h.

zlut_res

int CPUTracker::zlut_res

ZLUT resolution = number of radial steps in a plane.

Definition at line 64 of file cpu_tracker.h.

zluts

float* CPUTracker::zluts

Pointer to the first data in the ZLUTs.

All LUTs are saved in one big contiguous section of memory of size zlut_planes*zlut_count*zlut_res. Calculate specific LUTs or planes based on their indexes:

zluts[index * (zlut_planes * zlut_res) + plane * zlut_res + r]. 

The ZLUT system stores 'zlut_count' number of 2D zlut's, so every bead can be tracked with its own unique ZLUT.

Definition at line 61 of file cpu_tracker.h.

---

The documentation for this class was generated from the following files:

• cputrack/cpu_tracker.h
• cputrack/cpu_tracker.cpp