utils.h File Reference

#include "std_incl.h"
#include "scalar_types.h"

Go to the source code of this file.

Classes
struct	TImageData< T >
class	CImageData
struct	PathSeperator
class	Matrix3X3

Typedefs
typedef TImageData< float >	ImageData
typedef TImageData< double >	ImageDatad

Functions
template<typename T >
bool	isNAN (const T &v)
void	GetFormattedTimeString (char *output)
std::string	GetCurrentOutputPath (bool ext=true)
void	dbgout (const std::string &s)
std::string	SPrintf (const char *fmt,...)
void	dbgprintf (const char *fmt,...)
void	dbgsetlogfile (const char *path)
template<typename T >
void	DeleteAllElems (T &c)
template<typename TPixel >
void	normalize (TPixel *d, uint w, uint h)
template<typename T >
T	Lerp (T a, T b, float x)
template<typename T >
T	Interpolate (T *image, int width, int height, float x, float y, bool *outside=0)
template<typename T >
T	Interpolate1D (T *d, int len, float x)
template<typename T >
T	Interpolate1D (const std::vector< T > &d, float x)
template<typename T >
T	StdDeviation (T *start, T *end)
std::vector< float >	ComputeRadialBinWindow (int rsteps)
	Calculate the radial weights for ZLUT profile comparisons. More...
float	ComputeBgCorrectedCOM1D (float *data, int len, float cf=2.0f)
void	ComputeCRP (float *dst, int radialSteps, int angularSteps, float minradius, float maxradius, vector2f center, ImageData *src, float mean, float *crpmap=0)
void	ComputeRadialProfile (float *dst, int radialSteps, int angularSteps, float minradius, float maxradius, vector2f center, ImageData *src, float mean, bool normalize)
void	NormalizeRadialProfile (float *prof, int rsteps)
void	NormalizeZLUT (float *zlut, int numLUTs, int planes, int radialsteps)
void	GenerateImageFromLUT (ImageData *image, ImageData *zlut, float minradius, float maxradius, vector3f pos, bool useSplineInterp=true, int ovs=4)
void	ApplyPoissonNoise (ImageData &img, float poissonMax, float maxValue=255)
void	ApplyGaussianNoise (ImageData &img, float sigma)
void	WriteComplexImageAsCSV (const char *file, std::complex< float > *d, int w, int h, const char *labels[]=0)
void	WriteArrayAsCSVRow (const char *file, float *d, int len, bool append)
void	WriteVectorAsCSVRow (const char *file, std::vector< float > d, bool append)
void	WriteImageAsCSV (const char *file, float *d, int w, int h, const char *labels[]=0)
std::vector< std::vector< float > >	ReadCSV (const char *filename, char sep='\t')
std::vector< vector3f >	ReadVector3CSV (const char *file, char sep='\t')
void	WriteTrace (std::string file, vector3f *results, int nResults)
void	GenerateTestImage (ImageData &img, float xp, float yp, float size, float MaxPhotons)
std::string	GetLocalModuleFilename ()
std::string	GetLocalModulePath ()
std::string	GetDirectoryFromPath (std::string fullpath)
std::string	file_ext (const char *f)
ImageData	ReadJPEGFile (const char *fn)
ImageData	ReadLUTFile (const char *lutfile)
int	ReadJPEGFile (uchar *srcbuf, int srclen, uchar **data, int *width, int *height)
void	WriteJPEGFile (uchar *data, int w, int h, const char *filename, int quality)
void	FloatToJPEGFile (const char *name, const float *d, int w, int h)
void	WriteJPEGFile (const char *name, const ImageData &img)
int	NearestPowerOf2 (int v)
int	NearestPowerOf3 (int v)
void	GenerateGaussianSpotImage (ImageData *img, vector2f pos, float sigma, float I0, float Ibg)
std::vector< uchar >	ReadToByteBuffer (const char *filename)
double	GetPreciseTime ()
template<typename T >
void	floatToNormalizedInt (T *dst, const float *src, uint w, uint h, T maxValue)
template<typename T >
T *	floatToNormalizedInt (const float *src, uint w, uint h, T maxValue)
template<typename T >
T	ComputeStdDev (T *data, int len)
template<typename T >
T	qselect (T *data, int start, int end, int k)
template<typename T >
T	erf (T x)

Typedef Documentation

ImageData

typedef TImageData<float> ImageData

Definition at line 151 of file utils.h.

ImageDatad

typedef TImageData<double> ImageDatad

Definition at line 152 of file utils.h.

Function Documentation

ApplyGaussianNoise()

void ApplyGaussianNoise	(	ImageData &	img,
		float	sigma
	)

Definition at line 454 of file utils.cpp.

{
    for (int k=0;k<img.numPixels();k++) {
        float v = img.data[k] + sigma * rand_normal<float>();
        if (v<0.0f) v= 0.0f;
        img.data[k]=v;
    }
}

ApplyPoissonNoise()

void ApplyPoissonNoise	(	ImageData &	img,
		float	poissonMax,
		float	maxValue = 255
	)

Definition at line 432 of file utils.cpp.

{
    /*auto f = [&] (int y) {
        for (int x=0;x<img.w;x++) {
            img.at(x,y) = rand_poisson<float>(factor*img.at(x,y));
        }
    };

    ThreadPool<int, std::function<void (int index)> > pool(f);

    for (int y=0;y<img.h;y++) {
        pool.AddWork(y);
    }
    pool.WaitUntilDone();*/

    float ratio = maxval / poissonMax;

    for (int x=0;x<img.numPixels();x++) {
        img[x] = (int)(rand_poisson<float>(poissonMax*img[x]) * ratio );
    }
}

ComputeBgCorrectedCOM1D()

float ComputeBgCorrectedCOM1D	(	float *	data,
		int	len,
		float	cf = 2.0f
	)

Definition at line 231 of file utils.cpp.

{
    float sum=0, sum2=0;
    float moment=0;

    for (int x=0;x<len;x++) {
        float v = data[x];
        sum += v;
        sum2 += v*v;
    }

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

    for(int x=0;x<len;x++)
    {
        float v = data[x];
        v = std::max(0.0f, fabs(v-mean)-cf*stdev);
        sum += v;
        moment += x*v;
    }
    return moment / (float)sum;
}

ComputeCRP()

void ComputeCRP	(	float *	dst,
		int	radialSteps,
		int	angularSteps,
		float	minradius,
		float	maxradius,
		vector2f	center,
		ImageData *	src,
		float	mean,
		float *	crpmap = 0
	)

Definition at line 181 of file utils.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) );
    }

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

    float* map = crpmap ? crpmap : (float*)ALLOCA(sizeof(float)*radialSteps*angularSteps);
    float* com = (float*)ALLOCA(sizeof(float)*angularSteps);

    float rstep = (maxradius-minradius) / radialSteps;
    float comsum = 0.0f;
    for (int a=0;a<angularSteps;a++) {
        float r = minradius;
        float sum = 0.0f, moment=0.0f;
        for (int i=0;i<radialSteps; i++) {
            float x = center.x + radialDirs[a].x * r;
            float y = center.y + radialDirs[a].y * r;
            float v = img->interpolate(x,y);
            r += rstep;
            map[a*radialSteps+i] = v;
            sum += v;
            moment += i*v;
        }
        com[a] = moment/sum;
        comsum += com[a];
    }
    float avgcom = comsum/angularSteps;
    float totalrmssum2 = 0.0f;
    for (int i=0;i<radialSteps; i++) {
        double sum = 0.0f;
        for (int a=0;a<angularSteps;a++) {
            float shift = com[a]-avgcom;
            sum += map[a*radialSteps+i];
        }
        dst[i] = sum/angularSteps-paddingValue;
        totalrmssum2 += dst[i]*dst[i];
    }
    double invTotalrms = 1.0f/sqrt(totalrmssum2/radialSteps);
    for (int i=0;i<radialSteps;i++) {
        dst[i] *= invTotalrms;
    }
}

ComputeRadialBinWindow()

std::vector<float> ComputeRadialBinWindow

(

int

rsteps

)

Calculate the radial weights for ZLUT profile comparisons.

Calculates a mask, a so-called 'Stetson window', to reduce the effect of the edges of the profile in the calculation of error curves.

Parameters

[in]

rsteps

The number of radial steps in the window.

Definition at line 706 of file utils.cpp.

{
    std::vector<float> wnd(rsteps);
    for (int i=0;i<rsteps;i++) {
        float x = i/(float)rsteps;
        float t2 = 0.05f; 
        float t1 = 0.01f;
        float fall = 1.0f-expf(-sq(1-x)/t2);
        float rise = 1.0f-expf(-sq(x)/t1);
        wnd[i] = sqrtf(fall*rise*x*2);
    }
    return wnd;
}

ComputeRadialProfile()

void ComputeRadialProfile	(	float *	dst,
		int	radialSteps,
		int	angularSteps,
		float	minradius,
		float	maxradius,
		vector2f	center,
		ImageData *	src,
		float	mean,
		bool	normalize
	)

Definition at line 298 of file utils.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));
    }

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

//  center.x += 0.5f;
//  center.y += 0.5f;

    bool trace=false;
    float rstep = (maxradius-minradius) / radialSteps;
    int totalsmp = 0;

    //FILE* f = fopen("D:\\TestImages\\test.csv","w");

    for (int i=0;i<radialSteps; i++) {
        double sum = 0.0f;

        int nsamples = 0;
        float r = minradius+rstep*i;
        for (int a=0;a<angularSteps;a++) {
            float x = center.x + radialDirs[a].x * r;
            float y = center.y + radialDirs[a].y * r;
            bool outside;
            float v = img->interpolate(x,y, &outside);
            if (!outside) {
                sum += v;
                nsamples++;
            }

            //fprintf(f,"%d\t%d\t%f\n",i,a,v);
        }

        if (trace) {
            dbgprintf("%f,[%d]; ", sum, nsamples);
        }

        dst[i] = nsamples > MIN_RADPROFILE_SMP_COUNT ? sum/nsamples : mean;
    }
    if(trace)
        dbgprintf("\n");

    //fclose(f);

    if (normalize) 
        NormalizeRadialProfile(dst, radialSteps);
}

ComputeStdDev()

template<typename T >

T ComputeStdDev	(	T *	data,
		int	len
	)

Definition at line 221 of file utils.h.

{
    T sum = 0.0, sum2=0.0;
    for (int a=0;a<len;a++) {
        sum+=data[a];
        sum2+=data[a]*data[a];
    }
    T mean = sum / len;
    return sqrt(sum2 / len- mean * mean);
}

dbgout()

void dbgout

(

const std::string &

s

)

Definition at line 143 of file utils.cpp.

                                {
    OutputDebugString(s.c_str());
    printf(s.c_str());
    WriteToLog(s.c_str());
}

dbgprintf()

void dbgprintf	(	const char *	fmt,
			...
	)

Definition at line 149 of file utils.cpp.

                                    {
    va_list ap;
    va_start(ap, fmt);

    char buf[512];
    VSNPRINTF(buf, sizeof(buf), fmt, ap);
    OutputDebugString(buf);
    fputs(buf,stdout);
    WriteToLog(buf);

    va_end(ap);
}

dbgsetlogfile()

void dbgsetlogfile

(

const char *

path

)

Definition at line 49 of file utils.cpp.

{
    logFilename = path;
}

DeleteAllElems()

template<typename T >

void DeleteAllElems

(

T &

c

)

Definition at line 19 of file utils.h.

                          {
    for(typename T::iterator i=c.begin();i!=c.end();++i)
        delete *i;
    c.clear();
}

erf()

template<typename T >

T erf

(

T

x

)

Definition at line 373 of file utils.h.

{
    // constants
    T a1 =  0.254829592f;
    T a2 = -0.284496736f;
    T a3 =  1.421413741f;
    T a4 = -1.453152027f;
    T a5 =  1.061405429f;
    T p  =  0.3275911f;

    // Save the sign of x
    int sign = 1;
    if (x < 0)
        sign = -1;
    x = fabs(x);

    // A&S formula 7.1.26
    T t = 1.0f/(1.0f + p*x);
    T y = 1.0f - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return sign*y;
}

file_ext()

std::string file_ext

(

const char *

f

)

Definition at line 121 of file utils.cpp.

                                 {
    int l=strlen(f)-1;
    while (l > 0) {
        if (f[l] == '.')
            return &f[l+1];
        l--;
    }
    return "";
}

FloatToJPEGFile()

void FloatToJPEGFile	(	const char *	name,
		const float *	d,
		int	w,
		int	h
	)

Definition at line 189 of file fastjpg.cpp.

{
    uchar* zlut_bytes = floatToNormalizedInt(d, w,h, (uchar)255);
    WriteJPEGFile(zlut_bytes, w, h, name, 99);
    delete[] zlut_bytes;
}

floatToNormalizedInt() [1/2]

template<typename T >

void floatToNormalizedInt	(	T *	dst,
		const float *	src,
		uint	w,
		uint	h,
		T	maxValue
	)

Definition at line 199 of file utils.h.

{
    float maxv = src[0];
    float minv = src[0];
    for (uint k=0;k<w*h;k++) {
        maxv = std::max(maxv, src[k]);
        minv = std::min(minv, src[k]);
    }
    for (uint k=0;k<w*h;k++)
        dst[k] = maxValue * (src[k]-minv) / (maxv-minv);
}

floatToNormalizedInt() [2/2]

template<typename T >

T* floatToNormalizedInt	(	const float *	src,
		uint	w,
		uint	h,
		T	maxValue
	)

Definition at line 213 of file utils.h.

{ 
    T* r = new T[w*h]; 
    floatToNormalizedInt(r,src,w,h, maxValue);
    return r; 
}

GenerateGaussianSpotImage()

void GenerateGaussianSpotImage	(	ImageData *	img,
		vector2f	pos,
		float	sigma,
		float	I0,
		float	Ibg
	)

Definition at line 421 of file utils.cpp.

{
    float edenom = 1/sqrt(2*sigma*sigma);
    for (int y=0;y<img->h;y++)
        for(int x=0;x<img->w;x++) {
            float DeltaX = 0.5f * erf( (x-pos.x + .5f) * edenom ) - 0.5f * erf((x-pos.x - .5f) * edenom);
            float DeltaY = 0.5f * erf( (y-pos.y + .5f) * edenom ) - 0.5f * erf((y-pos.y - .5f) * edenom);
            img->at(x,y) = Ibg + I0 * DeltaX * DeltaY;
        }
}

GenerateImageFromLUT()

void GenerateImageFromLUT	(	ImageData *	image,
		ImageData *	zlut,
		float	minradius,
		float	maxradius,
		vector3f	pos,
		bool	useSplineInterp = true,
		int	ovs = 4
	)

Definition at line 354 of file utils.cpp.

{
//  lut.w = radialcov * ( (image->w/2 * roicov ) - minradius );
//  lut.w = radialcov * ( maxradius - minradius );

    float radialcov = zlut->w / (maxradius-minradius);
    float* zinterp = (float*)ALLOCA(zlut->w * sizeof(float));

    if (splineInterp) {
        int iz = std::max(1, std::min(zlut->h-3, (int)pos.z));
        float weights[4];
        float fz = pos.z-iz;
        ComputeBSplineWeights(weights, fz);
        // Interpolate ZLUT using B-spline weights
        for (int r=0;r<zlut->w;r++) {
            float zlutv = 0;
            for (int i=0;i<4;i++)
                zlutv += weights[i] * zlut->at(r, i-1+iz);
            zinterp[r] = zlutv;
        }
    }
    else {
        // The two Z planes to interpolate between
        int iz = (int)pos.z;
        if (iz < 0) 
            zinterp = zlut->data;
        else if (iz>=zlut->h-1)
            zinterp = &zlut->data[ (zlut->h-1)*zlut->w ];
        else {
            float* zlut0 = &zlut->data [ (int)pos.z * zlut->w ]; 
            float* zlut1 = &zlut->data [ ((int)pos.z + 1) * zlut->w ];
            zinterp = (float*)ALLOCA(sizeof(float)*zlut->w);
            for (int r=0;r<zlut->w;r++) 
                zinterp[r] = Lerp(zlut0[r], zlut1[r], pos.z-iz);
        }
    }

    int oversampleWidth=oversampleSubdiv,oversampleHeight=oversampleSubdiv;
    float oxstep = 1.0f / oversampleWidth;
    float oystep = 1.0f / oversampleHeight;

    for (int y=0;y<image->h;y++)
        for (int x=0;x<image->w;x++) 
        {
            float s = 0.0f;

            for (int ox=0;ox<oversampleWidth;ox++)
                for (int oy=0;oy<oversampleHeight;oy++) {

                    float X = x+(ox+0.5f)*oxstep - pos.x - 0.5f;
                    float Y = y+(oy+0.5f)*oystep - pos.y - 0.5f;

                    float pixr = sqrtf(X*X+Y*Y);
                    float r = (pixr - minradius) * radialcov;

                    if (r > zlut->w-2)
                        r = zlut->w-2;
                    if (r < 0) r = 0;

                    int i=(int)r;
                    s += Lerp(zinterp[i], zinterp[i+1], r-i);
                }
        
            image->at(x,y) = s/(oversampleWidth*oversampleHeight); 
        }
}

GenerateTestImage()

void GenerateTestImage	(	ImageData &	img,
		float	xp,
		float	yp,
		float	size,
		float	MaxPhotons
	)

Definition at line 162 of file utils.cpp.

{
    float S = 1.0f/sqrt(size);
    for (int y=0;y<img.h;y++) {
        for (int x=0;x<img.w;x++) {
            float X = x - xp;// + 0.5;
            float Y = y - yp;// + 0.5;
            float r = sqrtf(X*X+Y*Y)+1;
            float v = sinf(r/(5*S)) * expf(-r*r*S*0.001f);
            img.at(x,y)=v;
        }
    }

    if (SNratio>0) {
        ApplyGaussianNoise(img, 1.0f/SNratio);
    }
    img.normalize();
}

GetCurrentOutputPath()

std::string GetCurrentOutputPath

(

bool

ext = true

)

Definition at line 19 of file utils.cpp.

{
    std::string base = "D:\\TestImages\\TestOutput\\";
    std::string search = base + "*";
    LPCSTR w_folder = _T(search.c_str());
    
    WIN32_FIND_DATA FindFileData;
    HANDLE hFind;
    
    hFind = FindFirstFile(w_folder,&FindFileData);
    std::string dirName;
    while(FindNextFile(hFind,&FindFileData))
        dirName = FindFileData.cFileName;
    if(ext)
        return SPrintf("%s%s\\%s",base.c_str(),dirName.c_str(),"ZLUTDiag\\");
    else
        return SPrintf("%s%s",base.c_str(),dirName.c_str());
}

GetDirectoryFromPath()

std::string GetDirectoryFromPath

(

std::string

fullpath

)

Definition at line 104 of file utils.cpp.

{
    for (int i=fullpath.size()-1;i>=0;i--) {
        if (fullpath[i] == '/' || fullpath[i] == '\\') {
            return fullpath.substr(0, i);
        }
    }
    return "";
}

GetFormattedTimeString()

void GetFormattedTimeString

(

char *

output

)

Definition at line 38 of file utils.cpp.

{
    time_t rawtime;
    struct tm * timeinfo;
    time(&rawtime);
    timeinfo = localtime(&rawtime);
    sprintf(output, "%02d%02d%02d-%02d%02d%02d",timeinfo->tm_year-100,timeinfo->tm_mon+1,timeinfo->tm_mday,timeinfo->tm_hour,timeinfo->tm_min,timeinfo->tm_sec);
}

GetLocalModuleFilename()

std::string GetLocalModuleFilename

(

)

Definition at line 54 of file utils.cpp.

{
#ifdef WIN32
    char path[256];
    HMODULE hm = NULL;

    GetModuleHandleExA(GET_MODULE_HANDLE_EX_FLAG_FROM_ADDRESS | GET_MODULE_HANDLE_EX_FLAG_UNCHANGED_REFCOUNT,
                            (LPCSTR) &GetLocalModuleFilename, &hm);

    GetModuleFileNameA(hm, path, sizeof(path));
    return path;
#else
    #error GetLocalModuleName() not implemented for this platform
#endif
}

GetLocalModulePath()

std::string GetLocalModulePath

(

)

Definition at line 114 of file utils.cpp.

{
    std::string dllpath = GetLocalModuleFilename();
    return GetDirectoryFromPath(dllpath);
}

GetPreciseTime()

double GetPreciseTime

(

)

Definition at line 669 of file utils.cpp.

{
    uint64_t freq, time;

    QueryPerformanceCounter((LARGE_INTEGER*)&time);
    QueryPerformanceFrequency((LARGE_INTEGER*)&freq);

    return (double)time / (double)freq;
}

Interpolate()

template<typename T >

T Interpolate ( T *  image, int  width, int  height, float  x, float  y, bool *  outside = 0  )

inline

Definition at line 43 of file utils.h.

{
    int rx=x, ry=y;
    if (rx<0 || ry <0 || rx+1 >= width || ry+1>=height) {
        if (outside) *outside=true;
        return 0.0f;
    }
    if (outside) *outside=false;

    T v00 = image[width*ry+rx];
    T v10 = image[width*ry+rx+1];
    T v01 = image[width*(ry+1)+rx];
    T v11 = image[width*(ry+1)+rx+1];

    T v0 = Lerp(v00, v10, x-rx);
    T v1 = Lerp(v01, v11, x-rx);

    return Lerp(v0, v1, y-ry);
}

Interpolate1D() [1/2]

template<typename T >

T Interpolate1D ( T *  d, int  len, float  x  )

inline

Definition at line 64 of file utils.h.

{
    int fx = (int)x;
    if (fx < 0) return d[0];
    if (fx >= len-1) return d[len-1];
    return (d[fx+1]-d[fx]) * (x-fx) + d[fx];
}

Interpolate1D() [2/2]

template<typename T >

T Interpolate1D ( const std::vector< T > &  d, float  x  )

inline

Definition at line 73 of file utils.h.

{
    return Interpolate1D(&d[0],d.size(),x);
}

isNAN()

template<typename T >

bool isNAN

(

const T &

v

)

Definition at line 6 of file utils.h.

                                            { 
    return !(v == v); 
}

Lerp()

template<typename T >

T Lerp ( T  a, T  b, float  x  )

inline

Definition at line 40 of file utils.h.

{ return a + (b-a)*x; }

NearestPowerOf2()

int NearestPowerOf2

(

int

v

)

Definition at line 679 of file utils.cpp.

{
    int r=1;
    while (r < v) 
        r *= 2;
    if ( fabsf(r-v) < fabsf(r/2-v) )
        return r;
    return r/2;
}

NearestPowerOf3()

int NearestPowerOf3

(

int

v

)

Definition at line 689 of file utils.cpp.

{
    int r=1;
    while (r < v) 
        r *= 3;
    if ( fabsf(r-v) < fabsf(r/3-v) )
        return r;
    return r/3;
}

normalize()

template<typename TPixel >

void normalize	(	TPixel *	d,
		uint	w,
		uint	h
	)

Definition at line 27 of file utils.h.

{
    TPixel maxv = d[0];
    TPixel minv = d[0];
    for (uint k=0;k<w*h;k++) {
        maxv = std::max(maxv, d[k]);
        minv = std::min(minv, d[k]);
    }
    for (uint k=0;k<w*h;k++)
        d[k]=(d[k]-minv)/(maxv-minv);
}

NormalizeRadialProfile()

void NormalizeRadialProfile	(	float *	prof,
		int	rsteps
	)

Definition at line 257 of file utils.cpp.

{
    double sum=0.0f;
    for (int i=0;i<rsteps;i++)
        sum += prof[i];

    float mean =sum/rsteps;
    double rmssum2 = 0.0;

    for (int i=0;i<rsteps;i++) {
        prof[i] -= mean;
        rmssum2 += prof[i]*prof[i];
    }
    double invTotalrms = 1.0f/sqrt(rmssum2/rsteps);
    for (int i=0;i<rsteps;i++)
        prof[i] *= invTotalrms;
        
/*
    scalar_t minVal = prof[0];
    for (int i=0;i<rsteps;i++)
        if(prof[i]<minVal) minVal=prof[i]; 

    float rms=0;
    for (int i=0;i<rsteps;i++) {
        prof[i]-=minVal;
        rms += prof[i]*prof[i];
    }
    rms=1.0f/sqrt(rms);
    for (int i=0;i<rsteps;i++) {
        prof[i]*=rms;
    }*/
}

NormalizeZLUT()

void NormalizeZLUT	(	float *	zlut,
		int	numLUTs,
		int	planes,
		int	radialsteps
	)

Definition at line 291 of file utils.cpp.

{
    for(int i=0;i<numBeads;i++)
        for (int j=0;j<planes;j++)
            NormalizeRadialProfile(&zlut[radialsteps*planes*i + radialsteps*j], radialsteps);
}

qselect()

template<typename T >

T qselect	(	T *	data,
		int	start,
		int	end,
		int	k
	)

Definition at line 233 of file utils.h.

{
    if (end-start==1)
        return data[start];

    // select one of the elements as pivot
    int p = 0;
    T value = data[p+start];
    // swap with last value
    std::swap(data[p+start], data[end-1]);

    // move all items < pivot to the left
    int nSmallerItems=0;
    for(int i=start;i<end-1;i++)
        if(data[i]<value) {
            std::swap(data[i], data[start+nSmallerItems]);
            nSmallerItems++;
        }
    // pivot is now at [# items < pivot]
    std::swap(data[start+nSmallerItems], data[end-1]);

    // we are trying to find the kth element
    // so if pivotpos == k, we found it
    // if k < pivotpos, we need to recurse left side
    // if k > pivotpos, we need to recurse right side
    int pivotpos = start+nSmallerItems;
    if (k == pivotpos)
        return data[k];
    else if (k < pivotpos)
        return qselect(data, start, pivotpos, k);
    else 
        return qselect(data, pivotpos+1, end, k);
}

ReadCSV()

std::vector< std::vector<float> > ReadCSV	(	const char *	filename,
		char	sep = '\t'
	)

Definition at line 463 of file utils.cpp.

{
    std::list< std::vector <float> > data;

    FILE *f=fopen(filename,"r");
    if (f) {
        std::string buf;
        std::vector<float> vals;
        while (!feof(f)) {
            char c=fgetc(f);
            if (c == sep || c=='\n') {
                vals.push_back( atof(buf.c_str()) );
                buf.clear();
            }
            if (c == '\n') {
                data.push_back( vals );
                vals.clear();
            }

            if (c != sep && c!= '\n')
                buf+=c;
        }
        fclose(f);
    }
    
    std::vector<std::vector<float> > r;
    r.reserve(data.size());
    r.insert(r.begin(), data.begin(), data.end());
    return r;
}

ReadJPEGFile() [1/2]

ImageData ReadJPEGFile

(

const char *

fn

)

Definition at line 626 of file utils.cpp.

{
    int w, h;
    uchar* imgdata;
    std::vector<uchar> jpgdata = ReadToByteBuffer(fn);
    ReadJPEGFile(&jpgdata[0], jpgdata.size(), &imgdata, &w,&h);

    float* fbuf = new float[w*h];
    for (int x=0;x<w*h;x++)
        fbuf[x] = imgdata[x]/255.0f;
    delete[] imgdata;

    return ImageData(fbuf,w,h);
}

ReadJPEGFile() [2/2]

int ReadJPEGFile	(	uchar *	srcbuf,
		int	srclen,
		uchar **	data,
		int *	width,
		int *	height
	)

Definition at line 12 of file fastjpg.cpp.

{
  struct jpeg_decompress_struct cinfo;

  JSAMPARRAY buffer;        /* Output row buffer */
  int row_stride;       /* physical row width in output buffer */
  my_error_mgr jerr;
  cinfo.err = jpeg_std_error(&jerr.pub);
  jpeg_create_decompress(&cinfo);

  j_mem_src(&cinfo, srcbuf, srclen);

  /* Step 3: read file parameters with jpeg_read_header() */
  jpeg_read_header(&cinfo, TRUE);

  jpeg_start_decompress(&cinfo);
  row_stride = cinfo.output_width * cinfo.output_components;
  /* Make a one-row-high sample array that will go away when done with image */
  buffer = (*cinfo.mem->alloc_sarray) ((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1);

  *width = cinfo.output_width;
  *height = cinfo.output_height;
  *data = new uchar[cinfo.output_width*cinfo.output_height];

//  ResizeLVArray2D(output, cinfo.output_height, cinfo.output_width);

  /* Step 6: while (scan lines remain to be read) */
  /*           jpeg_read_scanlines(...); */

  /* Here we use the library's state variable cinfo.output_scanline as the
   * loop counter, so that we don't have to keep track ourselves.
   */
  uchar* dst = *data;
  while (cinfo.output_scanline < cinfo.output_height) {
    /* jpeg_read_scanlines expects an array of pointers to scanlines.
     * Here the array is only one element long, but you could ask for
     * more than one scanline at a time if that's more convenient.
     */
    jpeg_read_scanlines(&cinfo, buffer, 1);
    /* Assume put_scanline_someplace wants a pointer and sample count. */

    unsigned char* src = buffer[0];
    if (cinfo.output_components == 1) {
        memcpy(dst, src, cinfo.output_width);
    } else {
        for (uint x=0;x<cinfo.output_width;x++)
            dst[x] = src[x * cinfo.output_components];
    }
    dst += cinfo.output_width;
  }

  /* Step 7: Finish decompression */
  jpeg_finish_decompress(&cinfo);
  /* We can ignore the return value since suspension is not possible
   * with the stdio data source.
   */

  /* Step 8: Release JPEG decompression object */

  /* This is an important step since it will release a good deal of memory. */
  jpeg_destroy_decompress(&cinfo);

  /* After finish_decompress, we can close the input file.
   * Here we postpone it until after no more JPEG errors are possible,
   * so as to simplify the setjmp error logic above.  (Actually, I don't
   * think that jpeg_destroy can do an error exit, but why assume anything...)
   */
//  fclose(infile);

  /* At this point you may want to check to see whether any corrupt-data
   * warnings occurred (test whether jerr.pub.num_warnings is nonzero).
   */

  return 1;
}

ReadLUTFile()

ImageData ReadLUTFile

(

const char *

lutfile

)

Definition at line 720 of file utils.cpp.

{
    PathSeperator sep(lutfile);
    if(sep.extension == "jpg") {
        return ReadJPEGFile(lutfile);
    }
    else {
        std::string fn = lutfile;
        fn = std::string(fn.begin(), fn.begin()+fn.find('#'));
        std::string num( ++( sep.extension.begin() + sep.extension.find('#') ), sep.extension.end());
        int lutIndex = atoi(num.c_str());

        int nbeads, nplanes, nsteps;
        FILE *f = fopen(fn.c_str(), "rb");

        if (!f)
            throw std::runtime_error("Can't open " + fn);

        fread(&nbeads, 4, 1, f);
        fread(&nplanes, 4, 1, f);
        fread(&nsteps, 4, 1, f);


        fseek(f, 12 + 4* (nsteps*nplanes * lutIndex), SEEK_SET);
        ImageData lut = ImageData::alloc(nsteps,nplanes);
        fread(lut.data, 4, nsteps*nplanes,f);
        fclose(f);
        lut.normalize();
        return lut;
    }
}

ReadToByteBuffer()

std::vector<uchar> ReadToByteBuffer

(

const char *

filename

)

Definition at line 608 of file utils.cpp.

{
    FILE *f = fopen(filename, "rb");

    if (!f)
        throw std::runtime_error(SPrintf("%s was not found", filename));

    fseek(f, 0, SEEK_END);
    int len = ftell(f);
    fseek(f, 0, SEEK_SET);

    std::vector<uchar> buf(len);
    fread(&buf[0], 1,len, f);

    fclose(f);
    return buf;
}

ReadVector3CSV()

std::vector<vector3f> ReadVector3CSV	(	const char *	file,
		char	sep = '\t'
	)

Definition at line 494 of file utils.cpp.

{
    auto data=ReadCSV(file ,sep);

    std::vector<vector3f> r(data.size());

    for (uint i=0;i<data.size();i++){
        r[i]=vector3f(data[i][0],data[i][1],data[i][2]);
    }
    return r;
}

SPrintf()

std::string SPrintf	(	const char *	fmt,
			...
	)

Definition at line 132 of file utils.cpp.

                                        {
    va_list ap;
    va_start(ap, fmt);

    char buf[512];
    VSNPRINTF(buf, sizeof(buf), fmt, ap);

    va_end(ap);
    return buf;
}

StdDeviation()

template<typename T >

T StdDeviation	(	T *	start,
		T *	end
	)

Definition at line 139 of file utils.h.

                                 {
    T sum=0,sum2=0;
    for (T* s = start; s!=end; ++s) {
        sum+=*s; sum2+=(*s)*(*s);
    }

    T invN = 1.0f/(end-start);
    T mean = sum * invN;
    return sqrt(sum2 * invN - mean * mean);
}

WriteArrayAsCSVRow()

void WriteArrayAsCSVRow	(	const char *	file,
		float *	d,
		int	len,
		bool	append
	)

Definition at line 537 of file utils.cpp.

{
    FILE *f = fopen(file, append?"a":"w");
    if(f) {
        for (int i=0;i<len;i++)
            fprintf(f, "%.7f\t", d[i]);

        fprintf(f, "\n");
        fclose(f);
    }
    else
        dbgprintf("WriteArrayAsCSVRow: Unable to open file %s\n", file);
}

WriteComplexImageAsCSV()

void WriteComplexImageAsCSV	(	const char *	file,
		std::complex< float > *	d,
		int	w,
		int	h,
		const char *	labels[] = 0
	)

Definition at line 579 of file utils.cpp.

{
    FILE* f = fopen(file, "w");

    if (!f) {
        dbgprintf("WriteComplexImageAsCSV: Unable to open file %s\n", file);
        return;
    }

    if (labels) {
        for (int i=0;i<w;i++) {
            fprintf(f, "%s;\t", labels[i]);
        }
        fputs("\n", f);
    }

    for (int y=0;y<h;y++) {
        for (int x=0;x<w;x++)
        {
            float i=d[y*w+x].imag();
            fprintf(f, "%f%+fi", d[y*w+x].real(), i);
            if(x<w-1) fputs("\t", f); 
        }
        fprintf(f, "\n");
    }

    fclose(f);
}

WriteImageAsCSV()

void WriteImageAsCSV	(	const char *	file,
		float *	d,
		int	w,
		int	h,
		const char *	labels[] = 0
	)

Definition at line 551 of file utils.cpp.

{
    FILE* f = fopen(file, "w");

    if (f) {

        if (labels) {
            for (int i=0;i<w;i++) {
                fprintf(f, "%s;\t", labels[i]);
            }
            fputs("\n", f);
        }

        for (int y=0;y<h;y++) {
            for (int x=0;x<w;x++)
            {
                fprintf(f, "%.10f", d[y*w+x]);
                if(x<w-1) fputs("\t", f); 
            }
            fprintf(f, "\n");
        }

        fclose(f);
    }
    else
        dbgprintf("WriteImageAsCSV: Unable to open file %s\n", file);
}

WriteJPEGFile() [1/2]

void WriteJPEGFile	(	uchar *	data,
		int	w,
		int	h,
		const char *	filename,
		int	quality
	)

Definition at line 89 of file fastjpg.cpp.

{
  /* This struct contains the JPEG compression parameters and pointers to
   * working space (which is allocated as needed by the JPEG library).
   * It is possible to have several such structures, representing multiple
   * compression/decompression processes, in existence at once.  We refer
   * to any one struct (and its associated working data) as a "JPEG object".
   */
  struct jpeg_compress_struct cinfo;
  /* This struct represents a JPEG error handler.  It is declared separately
   * because applications often want to supply a specialized error handler
   * (see the second half of this file for an example).  But here we just
   * take the easy way out and use the standard error handler, which will
   * print a message on stderr and call exit() if compression fails.
   * Note that this struct must live as long as the main JPEG parameter
   * struct, to avoid dangling-pointer problems.
   */
  struct jpeg_error_mgr jerr;
  /* More stuff */
  JSAMPROW row_pointer[1];  /* pointer to JSAMPLE row[s] */
  int row_stride;       /* physical row width in image buffer */

  /* Step 1: allocate and initialize JPEG compression object */

  /* We have to set up the error handler first, in case the initialization
   * step fails.  (Unlikely, but it could happen if you are out of memory.)
   * This routine fills in the contents of struct jerr, and returns jerr's
   * address which we place into the link field in cinfo.
   */
  cinfo.err = jpeg_std_error(&jerr);
  /* Now we can initialize the JPEG compression object. */
  jpeg_create_compress(&cinfo);

  /* Step 2: specify data destination (eg, a file) */
  /* Note: steps 2 and 3 can be done in either order. */
  
  uint len = w*h * 2;
  uchar* memBuf = new uchar[len];
  j_mem_dest(&cinfo, (void**)&memBuf, &len);
  //jpeg_stdio_dest(&cinfo, outfile);

  /* Step 3: set parameters for compression */

  /* First we supply a description of the input image.
   * Four fields of the cinfo struct must be filled in:
   */
  cinfo.image_width = w;    /* image width and height, in pixels */
  cinfo.image_height = h;
  cinfo.input_components = 1;       /* # of color components per pixel */
  cinfo.in_color_space = JCS_GRAYSCALE;     /* colorspace of input image */
  /* Now use the library's routine to set default compression parameters.
   * (You must set at least cinfo.in_color_space before calling this,
   * since the defaults depend on the source color space.)
   */
  jpeg_set_defaults(&cinfo);
  /* Now you can set any non-default parameters you wish to.
   * Here we just illustrate the use of quality (quantization table) scaling:
   */
  jpeg_set_quality(&cinfo, quality, TRUE /* limit to baseline-JPEG values */);

  /* Step 4: Start compressor */

  /* TRUE ensures that we will write a complete interchange-JPEG file.
   * Pass TRUE unless you are very sure of what you're doing.
   */
  jpeg_start_compress(&cinfo, TRUE);

  /* Step 5: while (scan lines remain to be written) */
  /*           jpeg_write_scanlines(...); */

  /* Here we use the library's state variable cinfo.next_scanline as the
   * loop counter, so that we don't have to keep track ourselves.
   * To keep things simple, we pass one scanline per call; you can pass
   * more if you wish, though.
   */
  row_stride = w;   /* JSAMPLEs per row in image_buffer */

  while (cinfo.next_scanline < cinfo.image_height) {
    /* jpeg_write_scanlines expects an array of pointers to scanlines.
     * Here the array is only one element long, but you could pass
     * more than one scanline at a time if that's more convenient.
     */
    row_pointer[0] = &data[cinfo.next_scanline * row_stride];
    (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
  }
  jpeg_finish_compress(&cinfo);
  jpeg_destroy_compress(&cinfo);


  FILE *f = fopen(filename, "wb");
  if (f) {
      fwrite(memBuf, 1, len, f);
      fclose(f);
  } else 
      dbgprintf("Failed to open file %s for writing\n", filename);

  delete[] memBuf;
}

WriteJPEGFile() [2/2]

void WriteJPEGFile ( const char *  name, const ImageData &  img  )

inline

Definition at line 190 of file utils.h.

{ FloatToJPEGFile(name, img.data, img.w,img.h); }

WriteTrace()

void WriteTrace	(	std::string	file,
		vector3f *	results,
		int	nResults
	)

Definition at line 507 of file utils.cpp.

{
    FILE *f = fopen(filename.c_str(), "w");

    if (!f) {
        throw std::runtime_error(SPrintf("Can't open %s", filename.c_str()));
    }

    for (int i=0;i<nResults;i++)
    {
        fprintf(f, "%.7f\t%.7f\t%.7f\n", results[i].x, results[i].y, results[i].z);
    }

    fclose(f);
}

WriteVectorAsCSVRow()

void WriteVectorAsCSVRow	(	const char *	file,
		std::vector< float >	d,
		bool	append
	)

Definition at line 523 of file utils.cpp.

{
    FILE *f = fopen(file, append?"a":"w");
    if(f) {
        for (uint i=0;i<d.size();i++)
            fprintf(f, "%1.7f\t", d[i]);

        fprintf(f, "\n");
        fclose(f);
    }
    else
        dbgprintf("WriteArrayAsCSVRow: Unable to open file %s\n", file);
}