The math and image library provides a number of functions that are common to many modules in the pipeline, but not part of either the standard C library or GSL. This is where most of the heavy-lifting math is done. Math and image operations should generally be defined here with a generic interface rather than within a module so that they can be easily reused by other modules.

The function declarations are made within the name space math_func, and are defined in the file MathFunc.h. Any code that wishes to uses these function must include the header file, and either add a "using namespace math_func;" to the top of any function that uses math library functions, or provide the full name resolution when making the call, ie math_func:getRand().

Constants

double SNIM_PI
The value 3.14159265358979

double SNIM_ROOT_PI
The value 1.77245385090552

Basic Math Utilities

int roundToInt(double x)
Rounds x to the nearest integer. Works with both positive and negative values of x.

double factorial(int n)
Returns n! as a double value.

double nchoosek(int n, int k)
The number of ways to choose k objects from an n object set, without considering order. In other words, n!/(k!*(n-k)!).

double herm(int n, double x)
The value of the n-order cartesian hermite polynomial evaluated at x.

double hermPolar(int k, int l, double r)
The value of the k,l polar hermite polynomial evaluated at r.

Statistics

double covariance(double *d1, double *d2, int n)
Returns the covariance between two arrays of measurements. d1 and d2 are the arrays of values, n is the number of measurements.

Random Numbers

void setSeed(unsigned int s)
Sets the seed of the math_func random number generator to the given value. Note that this globally resets this pseudo-random sequence.

double getRand()
Returns a random double in the interval [0.0 1.0], using the math_func random number generator.

int getRandInt()
Returns a integer value in which every bit is random. The range is approximately -2^31 to 2^31.

unsigned int genRandomSeed()
Generates a completely random seed suitable for seeding pseudo-random sequences. Use of this function is preferable to seeding a sequence using time(), as every subsequent call generates a new seed, whereas time() will return the same number if a second has not passed.

Matrix Operations

double sum(double *m, int n)
Returns to sum of all elements of an array. m is a pointer to the array, and n is the number of elements in the array.

double sumProd(double *m1, double *m2, int n)
Returns the sum of the element-by-element product of the arrays. m1 and m2 are pointers to the arrays, and n is the number of elements in each array (which must be equal).

Geometrical Transformations and Utilities

void centroids(double *im, int w, int h, int numIms, double *flux, double *cx, double *cy)
Finds the centroids of a number of images by calculating their 1st order moments. im is a pointer to the image stack, consisting of numIms images of size w*h. flux, cx and cy are arrays of at least length numIms that will be filled in with the results of the calculation. flux will hold the total flux (sum) of each image. cx and cy will hold the x and y coordinates of the centroids for each image.

void coordTransform(double *t, double *xin, double *yin, double *xout, double *yout, int numel)
Takes a list of coordinate pairs, and performs the coordinate transform provided in t to each pair. t is the transform matrix (a four-element array), given as [t[0] t[2] ; t[1] t[3]]. xin and yin are arrays of the values of the input coordinates. xout and yout are arrays of at least equal size where the results of the calculation will be stored. numel is the number of elements in these array. An "in-place" transformation can be used by calling the function with xout=xin and yout=yin -- this overwrites the input values with the output values, using no additional memory.

Basis Transformations

void shapeletBasis2D(double *im, int w, int h, int n1, int n2, double cX, double cY, double beta)
Rasterize a cartesian shapelet base onto the provided image by sampling the basis function once per pixel. im is a double array of dimensions w and h onto which the basis function will be draw. n1 and n2 specify the order of the shapelet base. cX and cY specify the center of shapelet basis, in units of pixels where the cX = cY = 0 is the center of the first pixel of the image. beta specifies the scale function in units of pixels.

void shapeletBasisPolar2D(double *real, double *imaginary, int w, int h, int nl, int nr, double cX, double cY, double beta)
Rasterize a polar shapelet base onto the provided image by sampling the basis function once per pixel. im is a double array of dimensions w and h onto which the basis function will be draw. nl and nr specify the order of the shapelet base. cX and cY specify the center of shapelet basis, in units of pixels where the cX = cY = 0 is the center of the first pixel of the image. beta specifies the scale function in units of pixels.

void shapeletDecomp(double *coef, double *im, int w, int h, int numCoef, double cX, double cY, double beta)
Decompose an image into their cartesian shapelet coefficients. coef is a pointer to a double array in which the calculated coefficients will be stored, and must be at least numCoef*numCoef elements large. im is a pointer to the image that will be decomposed, and has dimensions w and h. The image will be decomposed up to numCoef in both n1 and n2 orders. cX and cY specify the center about which to do the decomposition, in units of pixels off set from the center of the first pixel (0,0). beta specifies the scale function to use in units of pixels.

void shapeletDecompSubsample(double *coef, double *im, int w, int h, int numCoef, double cX, double cY, double beta, int numSamples)
Perform the same operation with the same parameters as shapeletDecomp, but subsample the basis functions while performing the decomposition. numSamples*numSamples subsamples will be used for each pixel. If numSamples=1, the output will be identical to shapeletDecomp.

void shapeletRecon(double *coef, double *im, int w, int h, int numCoef, double cX, double cY, double beta)
Reconstruct an image from a set of cartesian shapelet coefficients. im is a pointer to an image array that will be filled in with the reconstructed values, and must be of size w*h. coef is a pointer to the coefficients to use, and must be of size numCoef*numCoef. cX and cY specify the center about which to do the reconstruction, in pixel units off set from the center of the first pixel (0,0). beta specifies the scale function in units of pixels.

void shapeletReconSubsample(double *coef, double *im, int w, int h, int numCoef, double cX, double cY, double beta, int numSubsamples)
Perform the same operation with the same parameters as shapeletRecon, but subsample the basis function while performing the reconstruction. numSubsamples*numSubsamples subsamples will be used for each pixel. If numSubsamples=1, the output will be identical to shapeletRecon.

Shear Estimators

void shapeletShear(double *meanGalCoef, double *coef, int numCoef, int numDat, double &g1, double &g2)
Calculate the shear bias in a group of coefficent arrays compared to a given average coefficient array. Typically, the latter array is unbiased on average compared to the former group. meanGalCoef is an array of the average coefficients, which must be of size numCoef*numCoef. The coef array contains numDat number of numCoef*numCoef shapelet coefficients. The pass-by-reference doubles g1 and g2 will be set to the result of the shear measurement.

void ksbMoments(double *im, int w, int h, int numIm, double *cx, double *cy, double *Qxx, double *Qxy, double *Qyy, double *e1, double *e2, Weight *weight)

double ksbShear(double *e, int n)
Calculate an estimate for shear bias (gamma_i) from the KSB elliptical moments (e_i). Namely, the function returns gamma_i = sum(e)/n. e is an array of elliptical moment measurements, and has length n.

void rrgMoments(double *im, int w, int h, int numIm, double *cx, double *cy, double *e1, double *e2, double *mu1, double *mu2, double *lam, double omega, Weight *weight);

void rrgShear(double *e1, double *e2, double *mu1, double *mu2, double *lam, int n, double &g1, double &g2)

Cosmic Ray Utilities