MathLib routines

/*MathLib routines

* MathLib.c * Please note last saved date below:> * $Date: 2008/05/16 16:04:19 $

***************************************************************************** * * Mark S. Cohen, Ph.D. c1999-2002 * All Rights reserved, except as noted below * * Unlimited academic use of this software is granted, except that any * publications that derive from its use, or any new software that * incorporates its underlying sources or algorithms must specifically * acknowledge the author. * * Modification of the program to meet local needs is encouraged, but all * such changes are to be communicated to the author to support the continued * growth of these tools. * * Commercial, for-profit, use of the sources or algorithms is expressly * prohibited without permission of the author. * */ /*************************************************************************** * MathLib contains a variety of routines that are largely for the * processing of complex data. *********************************************************************** * Keep control of MathLib revisions * void GetVLibVer( char *buff ) * * Returns the version number for MathLib in buff * ************************************************************************/ /* cvmov_new * Complex vector move void cvmov_new(float *source, int s_stride, float *dest, int d_stride, long count); /* ******************************************************************** TwoDGauss * Create a 2D gaussian kernel as a matrix of floats * INPUTS * data storage for output * xwidth, ywidth half-width (in pixels) of gaussian for x and y * xs, ys Matrix size. data must be large enough! * OUTPUTS * errorStatus * CHANGES * contents of data OSErr TwoDGauss( float *data, float xwidth, float ywidth, short xs, short ys ) *********************************************************************** Complex2DGauss * Create a 2D gaussian kernel as a matrix of floats * INPUTS * data storage for output * xwidth, ywidth half-width (in pixels) of gaussian for x and y * xs, ys Matrix size. data must be large enough! * OUTPUTS * errorStatus * CHANGES * contents of data * OSErr Complex2DGauss( float *data, float xwidth, float ywidth, short xs, short ys ) *********************************************************************** Complex3DGauss * Create a 3D gaussian kernel as a matrix of floats * INPUTS * data - storage for output * xwidth, ywidth, zwidth - half-width (in pixels) of gaussian for x and y * xs, ys, zs - Matrix size. data must be large enough! * OUTPUTS * errorStatus * CHANGES * contents of data * OSErr Complex3DGauss( float *data, float xwidth, float ywidth, ************************************************************************ vconvolve1D Function: vconvolve1D( float *invec, int vlength, float *outvec, float *kernel, int kernelLength ) Description: forms the convolution of invec with kernel and places it into outvec Parameters: invec - input data (float only) vlength - number of points in input data outvec - storage for output kernel - convolution kernel kernelLength - points in convolution kernel *********************************************************************** Boolean fft_bitrev_sincos_lut_init Function: fft_bitrev_sincos_lut_init() Description: Generate the lookup tables for fft functions Parameters: IN/OUT/IO [] Globals Accessed: bitrev2, ..., bitrev1024, sc_lut2, ..., sc_lut1024 Returns: completion flag. Boolean fft_bitrev_sincos_lut_init() *********************************************************************** bitreverse Function: bitreverse(int nsize, COMPLEX_F *pptr_cf, int *ppbitrev) Description: performs the bitreverse of the input vector Parameters: IN/OUT/IO [] nsize - size of the input vector *pptr_cf - pointer of the input vector with the struct COMPLEX_F *ppbitrev - pointer of the lookup table for the bitreverse Returns: error status - currently unimplemented. OSErr bitreverse(int nsize, COMPLEX *pptr_cf, int *ppbitrev) *********************************************************************** bitrev_norm Function: bitrev_norm(int nsize, COMPLEX_F *pptr_cf, int *ppbitrev) Description: performs the bitreverse and the normalization of the vector Parameters: IN/OUT/IO [] nsize - size of the input vector *pptr_cf - pointer of the input vector with the struct COMPLEX_F *ppbitrev - pointer of the lookup table for the bitreverse Returns: error status - currently unimplemented. OSErr bitrev_norm(int nsize, COMPLEX *pptr_cf, int *ppbitrev) *********************************************************************** cfft Function: cfft(float *buffer, int size, int direction) Description: In-place complex fft. Parameters: IN/OUT/IO [] IN buffer - pointer to complex array; size - input element count. direction: 1 for forward FFT; -1 for inverse FFT. OUT Places FFT output on top of input COMPLEX array. Globals Accessed: local global vectors bitrev and sc_lut Returns: error status. OSErr cfft(float *data, int size, char direction) *********************************************************************** cfft2d Function: cfft2d(float *ptr_f, int nx, int ny, char direction) Description: In-place 2D complex fft. Parameters: IN ptr_f- pointer to complex array; nx - input columns of elements; ny - input rows of elements; direction - 1 for forward FFT; -1 for inverse FFT. OUT Places 2D FFT output on top of input COMPLEX array. OSErr cfft2d(float *data, int nx, int ny, char direction) *********************************************************************** cfft3d Function: cfft3d(float *ptr_f, int nx, int ny, char direction) Description: In-place 3D complex fft. Parameters: IN ptr_f- pointer to complex array; nx - input columns of elements; ny - input rows of elements; nz - input planes of elements; direction - 1 for forward FFT; -1 for inverse FFT. OUT Places 3D FFT output on top of input COMPLEX array. How? Treat as a series of one-D transforms. Copy each line somewhere, do an FT then copy it back. Not wildly efficient: copies the data back and forth prior to and after each fft. Therefore, it is not a truly "in-place" routine. Fix someday? OSErr cfft3d(float *data, int nx, int ny, int nz, char direction) *********************************************************************** exch_quad Function : exch_quad Description : This routine will exchange the first and third quadrants, and exchange the second and fourth quadrants. The exchanged data will be found in the same buffer as the input data. Parameters : input_data - data to be quadrant swapped ysize - number of points in a full column xsize - number of points in a full row Returns : None void exch_quad( float *inbuf, short xsize, short ysize) *********************************************************************** exch_quad3d Function : exch_quad3d Description : This routine will exchange the corners in a complex 3d Volume. Parameters : input_data - data to be quadrant swapped ysize - number of points in a full column xsize - number of points in a full row zsize - number of planes Returns : None OSErr exch_quad3d( float *inbuf, short xsize, short ysize, short zsize) *********************************************************************** exch_vec void exch_vec( float *inbuf, long len ); /* *********************************************************************** vreal and vimag Function : vreal and vimag Description : Extract real or imaginary part of an input vector Parameters : inbuf - complex input vector outbuf - storage for output npairs - number of complex pairs in the inputvector Returns : None OSErr vreal( float *inbuf, float *outbuf, long nPairs ) OSErr vimag( float *inbuf, float *outbuf, long nPairs ) *********************************************************************** RealToComplex * void RealToComplex( float *vec, long npts ); * Takes a vector of real values and expands it to a complex vector * INPUTS * vec the input vector - must be large enough to handle complex size * npts number of real points * * OUTPUTS * none void RealToComplex( float *vec, long npts ) *********************************************************************** * vmag * Extract the magnitude from a complex vector - ok for in-place use * INPUTS * in the input vector * inStride number of floats separating adjacent complex pairs * outStride number of floats by which to separate magnitudes * nPairs * OUTPUTS * out the output vector * * OSErr vmag( float *in, short inStride, float *out, short outStride, long nPairs ) *********************************************************************** vphase * OSErr vphase( float *in, short inStride, float *out, short outStride, long nPairs); * Extract the phase (in radians) from a complex vector - ok for in-place use * INPUTS * in the input vector * inStride number of floats separating adjacent complex pairs * outStride number of floats by which to separate phases * nPairs * OUTPUTS * out the output vector * OSErr vphase( float *in, short inStride, float *out, short outStride, long nPairs ) *********************************************************************** vdegree * OSErr vdegree( float *in, short inStride, float *out, short outStride, long nPairs ) * Convert a vector in radians to a vector in degrees - ok for in-place use * INPUTS * in the input vector * inStride number of floats separating adjacent phases in radians * outStride number of floats by which to separate phases in degrees * nPairs * OUTPUTS * out the output vector *********************************************************************** * PadVolume * * This function performs padding or unpadding of image data as an * in-place operation. * * Space allocated for the padded image data must be * larger than the image size, because Fourier transform requires * data to have power of 2 elements in each dimension. If PAD is * selected, this function will pad the data evenly with zeros in all * three dimensions; if UNPAD is selected, it will move the data back * up to the beginning of the dataset. * * IN: theVolume - pointer to the input image data * im - UC image type for image dimensions * volXS, volYS, volZS - size of data as processed - power of 2 * * direction - PAD(1) or UNPAD(0); * * OUT: outVolume - desired location of output data * * * |<------------- fullXS ------------>| * |<- osX ->| | * padded --->.....................................----------- * ..................................... ^ * ..................................... osY | * .....................................---- | * ..........Soooooooooooooo............ | * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ fullYS * ..........ooooooooooooooo............ * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ | * ..........ooooooooooooooo............ | * ..................................... | * ..................................... | * ..................................... | * ..................................... V * .....................................----------- * * S = paddedStart * * E.g., * inXS = inYS = inZS = 4; * fullXS = fullYS = fullZS = 8; * * 3D padding (2X each dimension) * * Input: * Sxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * * * Output: * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ..Sxxx.. ..xxxx.. ..xxxx.. ..xxxx.. ........ ........ * ........ ........ ..xxxx.. ..xxxx.. ..xxxx.. ..xxxx.. ........ ........ * ........ ........ ..xxxx.. ..xxxx.. ..xxxx.. ..xxxx.. ........ ........ * ........ ........ ..xxxx.. ..xxxx.. ..xxxx.. ..xxxx.. ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * ........ ........ ........ ........ ........ ........ ........ ........ * * *********************************************************************** correlation * OSErr correlation( float *X, float *Y, double *r, double *slope, long npairs ) * Calculates the correlation between vector X and Y. Returns correlation, significance * and slope,b: v1 = b*v2. * INPUTS * X, Y input vectors * npairs Number of points in X and Y * OUTPUTS * r correlation (r, not r^2) of X and Y * slope slope of the regression of X and onto Y *********************************************************************** pcorr * OSErr pcorr( double r, double slope, double *p, double *conf95, long npairs ) * Calculates the probability of non-zero correlation. * INPUTS * r, calculated correlation * npairs Number of support points for r * OUTPUTS * p probability that v1 and v2 are correlated significantly * conf95 95% confidence interval for slope (can be omitted) * Assumes large n ( > 100 )