Index: /trunk/stac/src/Makefile.median =================================================================== --- /trunk/stac/src/Makefile.median (revision 3613) +++ /trunk/stac/src/Makefile.median (revision 3613) @@ -0,0 +1,35 @@ +SHELL = /bin/sh +CC = gcc +CFLAGS += -g -pg -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE -DCRFLUX -DTESTING +PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lcfitsio -lm -lxml2 -lmysqlclient +LDFLAGS += $(PSLIB) + +TARGET = median + +.PHONY: tags clean empty test profile optimise + +.c.o: + $(CC) -c $(CFLAGS) $(OPTFLAGS) $< + +median: median.o + $(CC) $(CFLAGS) -o $@ median.o $(LDFLAGS) $(OPTFLAGS) + +clean: + -$(RM) *.o gmon.* profile.txt + +empty: clean + -$(RM) $(TARGET) TAGS + +test: median + ./median test_3.fits.shift.1 + +# Run profiling. +profile: CFLAGS += -pg +profile: empty $(TARGET) + $(RM) gmon.* + for ((i = 0; i < 20; i++)) \ + do \ + $(MAKE) -f Makefile.median test; \ + mv -f gmon.out gmon.$$i; \ + done + gprof -p -q -l -c $(TARGET) gmon.* > profile.txt Index: /trunk/stac/src/fitTrans.c =================================================================== --- /trunk/stac/src/fitTrans.c (revision 3613) +++ /trunk/stac/src/fitTrans.c (revision 3613) @@ -0,0 +1,183 @@ +#include +#include "pslib.h" +#include "fitTrans.h" + +int main(int argc, char **argv) +{ + // Set trace levels + psTraceSetLevel(".", 0); + psTraceSetLevel("fitTrans.config", 10); + psTraceSetLevel("fitTrans.read", 10); + psTraceSetLevel("stac.checkMemory", 10); + + // Set logging level + psLogSetLevel(9); + + // Parse command line + fitTransConfig *config = fitTransParseConfig(argc, argv); + + printf("Order: %d\n", config->order); + + // Read the input files + psVector *xReference = NULL; + psVector *yReference = NULL; + psVector *xInput = NULL; + psVector *yInput = NULL; + (void)fitTransRead(&xReference, &yReference, config->reference); + (void)fitTransRead(&xInput, &yInput, config->input); + int nStars = xReference->n; // Number of stars + if (nStars != xInput->n) { + fprintf(stderr, "Number of stars in reference (%d) and input lists (%d) do not match.\n", + xReference->n, xInput->n); + exit(EXIT_FAILURE); + } + + psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates + psPlane *outCoord = psAlloc(sizeof(psPlane)); // Output Coordinates + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics + psVector *mask = psVectorAlloc(nStars, PS_TYPE_U8); // Mask for stars + for (int s = 0; s < nStars; s++) { + mask->data.U8[s] = 0; + } + +#if 0 + for (int i = 0; i < nStars; i++) { + xReference->data.F32[i] = (xReference->data.F32[i] - 512.0)/512.0; + yReference->data.F32[i] = (yReference->data.F32[i] - 512.0)/512.0; + xInput->data.F32[i] = (xInput->data.F32[i] - 512.0)/512.0; + yInput->data.F32[i] = (yInput->data.F32[i] - 512.0)/512.0; + } +#endif + + // Solve for the transformation + //psPlaneTransform *transform = fitTransSolve(NULL, xReference, yReference, xInput, yInput, NULL, config); + psPlaneTransform *transform = fitTransSolve(NULL, xReference, yReference, xInput, yInput, mask, config); + + // Iterate to do rejection + bool keepGoing = true; + for (int i = 0; i < config->nIter && keepGoing; i++) { + + psVector *xDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in x + psVector *yDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in y + for (int s = 0; s < nStars; s++) { + inCoord->x = xInput->data.F32[s]; + inCoord->y = yInput->data.F32[s]; + (void)psPlaneTransformApply(outCoord, transform, inCoord); + xDiff->data.F32[s] = xReference->data.F32[s] - outCoord->x; + yDiff->data.F32[s] = yReference->data.F32[s] - outCoord->y; + //printf("%f %f\n", xDiff->data.F32[s], yDiff->data.F32[s]); + } + + (void)psVectorStats(stats, xDiff, mask, 0); + float xSD = stats->sampleStdev; // Standard deviation in x + float xMean = stats->sampleMean; + (void)psVectorStats(stats, yDiff, mask, 0); + float ySD = stats->sampleStdev; // Standard deviation in y + float yMean = stats->sampleMean; + printf("Mean is %f, %f\n", xMean, yMean); + printf("Standard deviation is %f, %f\n", xSD, ySD); + + // Do the clipping + int nRej = 0; // Number rejected + for (int s = 0; s < nStars; s++) { + if (!mask->data.U8[s] && ((fabsf(xDiff->data.F32[s]) > config->nSigma * xSD) || + (fabsf(yDiff->data.F32[s]) > config->nSigma * ySD))) { + mask->data.U8[s] = 1; + nRej++; + printf(" Clipping star %d\n", s); + } + } + printf("%d stars clipped on iteration %d.\n", nRej, i + 1); + psFree(xDiff); + psFree(yDiff); + + if (nRej == 0) { + keepGoing = false; // Nothing more to reject + } else { + fitTransSolve(transform, xReference, yReference, xInput, yInput, mask, config); + } + } + + psVector *xDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in x + psVector *yDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in y + for (int s = 0; s < nStars; s++) { + inCoord->x = xInput->data.F32[s]; + inCoord->y = yInput->data.F32[s]; + (void)psPlaneTransformApply(outCoord, transform, inCoord); + xDiff->data.F32[s] = xReference->data.F32[s] - outCoord->x; + yDiff->data.F32[s] =yReference->data.F32[s] - outCoord->y; + } + + (void)psVectorStats(stats, xDiff, mask, 0); + float xSD = stats->sampleStdev; // Standard deviation in x + float xMean = stats->sampleMean; // Mean of x + (void)psVectorStats(stats, yDiff, mask, 0); + float ySD = stats->sampleStdev; // Standard deviation in y + float yMean = stats->sampleMean; // Mean of y + printf("Final mean is %f, %f\n", xMean, yMean); + printf("Final standard deviation is %f, %f\n", xSD, ySD); + psFree(xDiff); + psFree(yDiff); + + psFree(inCoord); + psFree(outCoord); + psFree(mask); + psFree(stats); + + int order = config->order + 1; // Polynomial order + int nCoeff = order * (order + 1) / 2; // Number of coefficients + printf("x' = "); + for (int i = 0; i < order; i++) { + for (int j = 0; j < order - i; j++) { + if (i == 0 && j == 0) { + printf ("%f\n", transform->x->coeff[0][0]); + } else { + printf (" + %f x^%d y^%d\n", transform->x->coeff[i][j], i, j); + } + } + } + printf("\n"); + printf("y' = "); + for (int i = 0; i < order; i++) { + for (int j = 0; j < order - i; j++) { + if (i == 0 && j == 0) { + printf ("%f\n", transform->y->coeff[0][0]); + } else { + printf (" + %f x^%d y^%d\n", transform->y->coeff[i][j], i, j); + } + } + } + + if (config->mapFile != NULL) { + FILE *fp = fopen(config->mapFile, "w"); // File pointer + if (fp == NULL) { + fprintf(stderr, "Unable to open %s for writing.\n", fp); + exit(EXIT_FAILURE); + } else { + fprintf(fp, "%d\n", config->order); + for (int i = 0; i < order; i++) { + for (int j = 0; j < order - i; j++) { + fprintf(fp, "%.10f ", transform->x->coeff[i][j]); + } + } + fprintf(fp, "\n"); + for (int i = 0; i < order; i++) { + for (int j = 0; j < order - i; j++) { + fprintf(fp, "%.10f ", transform->y->coeff[i][j]); + } + } + fprintf(fp, "\n"); + fclose(fp); + } + } + + // Clean up + psFree(xReference); + psFree(yReference); + psFree(xInput); + psFree(yInput); + psFree(transform); + psFree(config); + + stacCheckMemory(); +} Index: /trunk/stac/src/fitTrans.h =================================================================== --- /trunk/stac/src/fitTrans.h (revision 3613) +++ /trunk/stac/src/fitTrans.h (revision 3613) @@ -0,0 +1,77 @@ +#ifndef FIT_TRANS_H +#define FIT_TRANS_H + + +// Temporary implementation +psPlaneTransform *psPlaneTransformAlloc(psS32 n1, psS32 n2); + +/************************************************************************************************************/ + +// Configuration for the program +typedef struct { + int verbose; // Verbosity level + int order; // Polynomial order of the fit + float nSigma; // Number of standard deviations at which to clip + int nIter; // Number of rejection iterations + const char *reference; // Name of reference star list + const char *input; // Name of input star list + const char *mapFile; // Name of output map file +} fitTransConfig; + + +/************************************************************************************************************/ + +// fitTransConfig.c + +// The help message +void help(const char *name + ); + +// Parse the command line arguments +fitTransConfig *fitTransParseConfig(int argc, // Number of arguments + char **argv // Arguments + ); + +/************************************************************************************************************/ + +// fitTransRead.c + +// Read an input file +bool fitTransRead(psVector **x, // x positions + psVector **y, // y positions + const char *filename // Filename from which to read + ); + +/************************************************************************************************************/ + +// stacCheckMemory.c +// This stuff was written for STAC, but it should work anywhere, so I'll just plug it in here. + +// Check memory +void stacCheckMemory(void); + +// Print out the problem +void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block. + const char *file, ///< the file in which the problem originated + psS32 lineno ///< the line number in which the problem originated + ); + +// Print out a memblock when it's allocated --- this function used as a callback +psMemoryId stacMemPrint(const psMemBlock *ptr); + +/************************************************************************************************************/ + +// fitTransSolve.c + +// Solve for the transformation +psPlaneTransform *fitTransSolve(psPlaneTransform *trans,// Output transformation, or NULL + psVector *xReference, psVector *yReference, // Reference coordinates + psVector *xInput, psVector *yInput, // Input coordinates + psVector *mask, // Mask for coordinates + fitTransConfig *config // Configuration + ); + +/************************************************************************************************************/ + + +#endif Index: /trunk/stac/src/fitTransConfig.c =================================================================== --- /trunk/stac/src/fitTransConfig.c (revision 3613) +++ /trunk/stac/src/fitTransConfig.c (revision 3613) @@ -0,0 +1,92 @@ +#include +#include +#include +#include "pslib.h" +#include "fitTrans.h" + + +void help(const char *name // Program name + ) +{ + fprintf (stderr, "fitTrans: Fit image transformation given two matched star lists.\n" + "Usage: %s [-h] [-v] [-o ORDER] [-s SIGMA] [-n NITER] [-m MAPFILE] REFERENCE INPUT\n" + "where\n" + "\t-h Help (this info)\n" + "\t-v Increase verbosity level\n" + "\t-o ORDER Order for polynomial fit: 1, 2, 3 (1)\n" + "\t-s SIGMA Number of standard deviations at which to reject (3)\n" + "\t-n NITER Number of rejection iterations (1)\n" + "\t-m MAPFILE Name of file to which to write the transformation map\n" + "\tREFERENCE Reference list: index, x, y, magnitude\n" + "\tINPUT Input list: index, x, y, magnitude\n", + name + ); +} + + +fitTransConfig *fitTransParseConfig(int argc, // Number of arguments + char **argv // Arguments + ) +{ + const char *programName = argv[0]; // Program name + + // Variables for getopt + int opt; // Option, from getopt + extern char *optarg; // Argument accompanying switch + extern int optind; // getopt variable + + // Initialise the configuration + fitTransConfig *config = psAlloc(sizeof(fitTransConfig)); + config->verbose = 0; // Verbosity level + config->order = 1; // Order for fit + config->nSigma = 3.0; // Number of sigma at which to clip + config->nIter = 1; // Number of rejection iterations + config->mapFile = NULL; // Name of map file + + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvo:s:n:m:")) != -1) { + switch (opt) { + case 'h': + help(programName); + exit(EXIT_SUCCESS); + case 'v': + config->verbose++; + break; + case 'o': + if (sscanf(optarg, "%d", &config->order) != 1) { + help(programName); + } + break; + case 's': + if (sscanf(optarg, "%f", &config->nSigma) != 1) { + help(programName); + } + break; + case 'n': + if (sscanf(optarg, "%d", &config->nIter) != 1) { + help(programName); + } + break; + case 'm': + config->mapFile = optarg; + break; + default: + help(programName); + } + } + + argc -= optind; + argv += optind; + if (argc != 2) { + help(programName); + exit(EXIT_FAILURE); + } + + config->reference = argv[0]; // Reference list + config->input = argv[1]; // Input list + + psTrace("fitTrans.config", 1, "Reference file is %s\n", config->reference); + psTrace("fitTrans.config", 1, "Input file is %s\n", config->input); + + return config; +} Index: /trunk/stac/src/fitTransRead.c =================================================================== --- /trunk/stac/src/fitTransRead.c (revision 3613) +++ /trunk/stac/src/fitTransRead.c (revision 3613) @@ -0,0 +1,41 @@ +#include +#include "pslib.h" +#include "fitTrans.h" + +#define BUFFER_SIZE 100 + +bool fitTransRead(psVector **x, // x positions + psVector **y, // y positions + const char *filename // Filename from which to read + ) +{ + FILE *fp; // File pointer + + fp = fopen(filename,"r"); + if (fp == NULL) { + psError("fitTrans.read", "Unable to open file %s\n", filename); + return false; + } + + int index = 0; // Index of star + float magnitude = 0.0; // Magnitude of star + int num = 0; // Number of stars read + + *x = psVectorAlloc(BUFFER_SIZE, PS_TYPE_F32); + *y = psVectorAlloc(BUFFER_SIZE, PS_TYPE_F32); + + while (fscanf(fp, "%d %f %f %f", &index, &((*x)->data.F32[num]), &((*y)->data.F32[num]), &magnitude) != EOF) { + num++; + if (num % BUFFER_SIZE) { + (void)psVectorRealloc(*x, num+BUFFER_SIZE); + (void)psVectorRealloc(*y, num+BUFFER_SIZE); + } + } + + (*x)->n = num; + (*y)->n = num; + + psTrace("fitTrans.read", 1, "Read %d entries from %s\n", num, filename); + + return true; +} Index: /trunk/stac/src/fitTransSolve.c =================================================================== --- /trunk/stac/src/fitTransSolve.c (revision 3613) +++ /trunk/stac/src/fitTransSolve.c (revision 3613) @@ -0,0 +1,127 @@ +#include +#include +#include "pslib.h" +#include "fitTrans.h" + + +psPlaneTransform *fitTransSolve(psPlaneTransform *trans,// Output coefficients, or NULL + psVector *xReference, psVector *yReference, // Reference coordinates + psVector *xInput, psVector *yInput, // Input coordinates + psVector *mask, // Mask for coordinates + fitTransConfig *config // Configuration + ) +{ + int order = config->order + 1; // Order of polynomials + int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients + + assert(!trans || ((trans->x->nX == order) && (trans->x->nY == order) && (trans->y->nX == order) && + (trans->y->nY == order))); + assert(xReference); + assert(yReference); + assert(xReference->n == yReference->n); + assert(xInput); + assert(yInput); + assert(xInput->n == yInput->n); + assert(xReference->n == xInput->n); + assert(!mask || xReference->n == mask->n); + + // Check inputs + if (trans == NULL) { + trans = psPlaneTransformAlloc(order, order); + // Initialise transformation + for (int i = 0; i < order; i++) { + for (int j = 0; j < order; j++) { + trans->x->coeff[i][j] = 0.0; + trans->y->coeff[i][j] = 0.0; + if (i + j < order) { + trans->x->mask[i][j] = 0; + trans->y->mask[i][j] = 0; + } else { + trans->x->mask[i][j] = 1; + trans->y->mask[i][j] = 1; + } + } + } + } + + // Initialise calculations + psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution + psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x + psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y + for (int i = 0; i < nCoeff; i++) { + for (int j = 0; j < nCoeff; j++) { + matrix->data.F64[i][j] = 0.0; + } + xVector->data.F64[i] = 0.0; + yVector->data.F64[i] = 0.0; + } + + // Create fake polynomial to use in evaluation + psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD); + for (int i = 0; i < order; i++) { + for (int j = 0; j < order; j++) { + fakePoly->coeff[i][j] = 1.0;// Set all coeffecients to 1 + fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate + } + } + + // Iterate over stars that haven't been masked + int numStars = 0; // Number of stars used + for (int k = 0; k < xReference->n; k++) { + if (!mask || ! mask->data.U8[k]) { + numStars++; + + // Iterate over the polynomial coefficients + for (int i = 0, ijIndex = 0; i < order; i++) { + for (int j = 0; j < order - i; j++, ijIndex++) { + + fakePoly->mask[i][j] = 0; + double ijPoly = psDPolynomial2DEval((double)xInput->data.F32[k], (double)yInput->data.F32[k], fakePoly); + fakePoly->mask[i][j] = 1; + + for (int m = 0, mnIndex = 0; m < order; m++) { + for (int n = 0; n < order - m; n++, mnIndex++) { + + fakePoly->mask[m][n] = 0; + double mnPoly = psDPolynomial2DEval((double)xInput->data.F32[k], (double)yInput->data.F32[k], fakePoly); + fakePoly->mask[m][n] = 1; + + matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly; + } + } + + xVector->data.F64[ijIndex] += ijPoly * (double)xReference->data.F32[k]; + yVector->data.F64[ijIndex] += ijPoly * (double)yReference->data.F32[k]; + } + } // Iterating over coefficients + } + } // Iterating over stars + + printf("Solving based on %d stars...\n", numStars); + + // Solution via LU Decomposition + psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition + psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix); // LU decomposed matrix + psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x + psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y + psFree(permutation); + psFree(luMatrix); + + // Stuff coefficients into transformation + for (int i = 0, ijIndex = 0; i < order; i++) { + for (int j = 0; j < order - i; j++, ijIndex++) { + trans->x->coeff[i][j] = xSolution->data.F64[ijIndex]; + trans->y->coeff[i][j] = ySolution->data.F64[ijIndex]; + } + } + + // Clean up + psFree(xSolution); + psFree(ySolution); + psFree(matrix); + psFree(xVector); + psFree(yVector); + psFree(fakePoly); + + return trans; +} Index: /trunk/stac/src/median.c =================================================================== --- /trunk/stac/src/median.c (revision 3613) +++ /trunk/stac/src/median.c (revision 3613) @@ -0,0 +1,59 @@ +#include +#include "pslib.h" + +double getTime(void) +/* Gets the current time. Got this from Nick Kaiser's fetchpix.c */ +{ + struct timeval tv; + gettimeofday(&tv, NULL); + return(tv.tv_sec + 1.e-6 * tv.tv_usec); +} + +int main(int argc, char *argv[]) +{ +#if 0 + psTraceSetLevel(".", 10); + psTraceSetLevel(".psLib", 10); + psTraceSetLevel(".psLib.dataManip", 10); + psTraceSetLevel(".psLib.dataManip.psFunctions", 10); + psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0); + psTraceSetLevel(".psLib.dataManip.psFunctions.vectorBinDisect", 0); + psTraceSetLevel(".pslib.stats.robust", 10); +#endif + + if (argc != 2) { + exit(EXIT_FAILURE); + } + + const char *filename = argv[1]; + + psFits *imageFile = psFitsAlloc(filename); + psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read + psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0); + if (image == NULL) { + psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", filename); + exit(EXIT_FAILURE); + } + psTrace("median", 4, "Image %s is %dx%d\n", filename, image->numCols, image->numRows); + psFree(imageRegion); + // Convert to 32-bit floating point, in necessary + if (image->type.type != PS_TYPE_F32) { + psTrace("stac.read.images", 3, "Converting %s to floating point in memory....", filename); + psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32); + psFree(image); + image = temp; + } + psFree(imageFile); + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + double startTime = getTime(); + stats = psImageStats(stats, image, NULL, 0); + double endTime = getTime(); + printf("Median: %f\n", stats->robustMedian); + printf("Time: %f\n", endTime - startTime); + + psFree(stats); + psFree(image); + + exit(EXIT_SUCCESS); +} Index: /trunk/stac/src/photScale.c =================================================================== --- /trunk/stac/src/photScale.c (revision 3613) +++ /trunk/stac/src/photScale.c (revision 3613) @@ -0,0 +1,341 @@ +#include +#include +#include +#include "pslib.h" +#include "psMinimize.h" + +#define BUFFER 100 // Buffer size for star list +#define MAXITER 100 // Maximum iterations for fitting +#define TOLERANCE 1e-2 // Fitting tolerance +#define ZP 25 // Zero point magnitude + +#define SQUARE(x) ((x)*(x)) +#define MAX(x,y) ((x) > (y) ? (x) : (y)) +#define MIN(x,y) ((x) < (y) ? (x) : (y)) + + +psVector *gaussian2d_deriv(psImage *deriv, // Derivatives wrt the parameters + const psVector *params,// Parameters with which to evaluate + const psArray *position // Position at which to evaluate + ) +{ + float bg = params->data.F32[0]; + float norm = params->data.F32[1]; + float x0 = params->data.F32[2]; + float y0 = params->data.F32[3]; + float sigmaU = params->data.F32[4]; + float sigmaV = params->data.F32[5]; + float theta = params->data.F32[6]; + + psVector *gauss = psVectorAlloc(position->n, PS_TYPE_F32); + for (int i = 0; i < position->n; i++) { + psVector *pos = position->data[i]; + + float x = pos->data.F32[0]; + float y = pos->data.F32[1]; + + float xSub = x - x0; + float ySub = y - y0; + + float cosTheta = cos(theta); + float sinTheta = sin(theta); + + float invSigmaU2 = 1.0 / SQUARE(sigmaU); + float invSigmaV2 = 1.0 / SQUARE(sigmaV); + + float u = - xSub * cosTheta + ySub * sinTheta; + float v = xSub * cosTheta + ySub * sinTheta; + + float uOnSigmaU2 = u * invSigmaU2; + float vOnSigmaV2 = v * invSigmaV2; + + float exponent = SQUARE(u) * invSigmaU2 + SQUARE(v) * invSigmaV2; + float bell = expf(- 0.5 * exponent); + gauss->data.F32[i] = norm * bell + bg; + + float deriv_bg = 1.0; + float deriv_norm = bell; + float deriv_x0 = norm * bell * (- uOnSigmaU2 + vOnSigmaV2) * cosTheta; + float deriv_y0 = norm * bell * (uOnSigmaU2 + vOnSigmaV2) * sinTheta; + float deriv_sigmaU = norm * bell * uOnSigmaU2 * u / sigmaU; + float deriv_sigmaV = norm * bell * vOnSigmaV2 * v / sigmaV; + float deriv_theta = - norm * bell * (uOnSigmaU2 * (xSub * sinTheta + ySub * cosTheta) + + vOnSigmaV2 * (-xSub * sinTheta + ySub * cosTheta)); + + deriv->data.F32[i][0] = deriv_bg; + deriv->data.F32[i][1] = deriv_norm; + deriv->data.F32[i][2] = deriv_x0; + deriv->data.F32[i][3] = deriv_y0; + deriv->data.F32[i][4] = deriv_sigmaU; + deriv->data.F32[i][5] = deriv_sigmaV; + deriv->data.F32[i][6] = deriv_theta; + } + + return gauss; +} + +psVector *gaussian2d_noderiv(const psVector *params,// Parameters with which to evaluate + const psArray *position // Position at which to evaluate + ) +{ + float bg = params->data.F32[0]; + float norm = params->data.F32[1]; + float x0 = params->data.F32[2]; + float y0 = params->data.F32[3]; + float sigmaU = params->data.F32[4]; + float sigmaV = params->data.F32[5]; + float theta = params->data.F32[6]; + + psVector *gauss = psVectorAlloc(position->n, PS_TYPE_F32); + for (int i = 0; i < position->n; i++) { + psVector *pos = position->data[i]; + + float x = pos->data.F32[0]; + float y = pos->data.F32[1]; + + float xSub = x - x0; + float ySub = y - y0; + + float cosTheta = cos(theta); + float sinTheta = sin(theta); + + float invSigmaU2 = 1.0 / SQUARE(sigmaU); + float invSigmaV2 = 1.0 / SQUARE(sigmaV); + + float u = - xSub * cosTheta + ySub * sinTheta; + float v = xSub * cosTheta + ySub * sinTheta; + + float uOnSigmaU2 = u * invSigmaU2; + float vOnSigmaV2 = v * invSigmaV2; + + float exponent = SQUARE(u) * invSigmaU2 + SQUARE(v) * invSigmaV2; + float bell = expf(- 0.5 * exponent); + gauss->data.F32[i] = norm * bell + bg; + + } + + return gauss; +} + + + +void help(char *programName) +{ + fprintf(stderr, "photScale: calculates the photometric scaling between two images\n" + "Usage: %s [-r REJECT] [-i NUM_ITER] REFERENCE INPUT STARS\n" + "where\n" + "\t-r REJECT Rejection limit, standard deviations (2.5)\n" + "\t-i NUM_ITER Number of rejection iterations (3)\n" + "\tREFERENCE Reference image\n" + "\tINPUT Input image\n" + "\tSTARS File with stellar positions, x y on each line\n", + programName); +} + + +int main(int argc, char *argv[]) +{ + int verbose = 0; // Verbosity level + int numIter = 3; // Number of iterations + float reject = 2.5; // Rejection limit (sigma) + int size = 5; // Size of stamps + + psTraceSetLevel(".", 10); + psTraceSetLevel("read", 10); + psTraceSetLevel("phot", 10); + psTraceSetLevel("gauss", 9); + psTraceSetLevel(".psLib.dataManip.psMinimize", 0); + + + /* Variables for getopt */ + int opt; /* Option, from getopt */ + extern char *optarg; /* Argument accompanying switch */ + extern int optind; /* getopt variables */ + + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvr:i:s:")) != -1) { + switch (opt) { + case 'h': + help(argv[-1]); + exit(EXIT_SUCCESS); + case 'v': + verbose++; + break; + case 'r': + if (sscanf(optarg, "%f", &reject) != 1) { + help(argv[-1]); + exit(EXIT_FAILURE); + } + break; + case 'i': + if (sscanf(optarg, "%d", &numIter) != 1) { + help(argv[-1]); + exit(EXIT_FAILURE); + } + break; + case 's': + if (sscanf(optarg, "%d", &size) != 1) { + help(argv[-1]); + exit(EXIT_FAILURE); + } + default: + help(argv[-1]); + exit(EXIT_FAILURE); + } + } + + /* Get the remaining, mandatory values */ + argc -= optind; + argv += optind; + if (argc != 3) { + help(argv[-1]); + exit(EXIT_FAILURE); + } + + const char *reference = argv[0]; + const char *input = argv[1]; + const char *stars = argv[2]; + + // Read images + psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read + psFits *referenceFile = psFitsAlloc(reference); // FITS file + psImage *refImage = psFitsReadImage(NULL, referenceFile, *imageRegion, 0); // Reference image + if (refImage == NULL) { + psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", reference); + exit(EXIT_FAILURE); + } + if (refImage->type.type != PS_TYPE_F32) { + // Convert to F32 + psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32); + psFree(refImage); + refImage = temp; + } + psFree(referenceFile); + psImageClipNaN(refImage, 0.0); + +#if 0 + psFits *inputFile = psFitsAlloc(input); // FITS file + psImage *inImage = psFitsReadImage(NULL, inputFile, *imageRegion, 0); // Input image + if (inImage == NULL) { + psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", input); + exit(EXIT_FAILURE); + } + if (inImage->type.type != PS_TYPE_F32) { + // Convert to F32 + psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32); + psFree(inImage); + refImage = temp; + } + psFree(inputFile); +#endif + + // List of star positions + psArray *starList = psArrayAlloc(BUFFER); + FILE *starFile = fopen(stars, "r"); + if (starFile == NULL) { + fprintf(stderr, "Unable to open %s\n", stars); + exit(EXIT_FAILURE); + } + float x, y; + int num = 0; // Number of stars + while (fscanf(starFile, "%f %f", &x, &y) == 2) { + psPlane *coords = psPlaneAlloc(); + coords->x = x; + coords->y = y; + starList->data[num] = coords; + num++; + if (num % BUFFER) { + starList = psArrayRealloc(starList, num + BUFFER); + } + } + starList->n = num; + psTrace("read", 1, "%d stars read.\n", num); + + psStats *medianStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + (void)psImageStats(medianStats, refImage, NULL, 0); + float refMedian = medianStats->sampleMedian; + + psTrace("phot", 3, "Image median: %f\n", refMedian); + + psVector *refFlux = psVectorAlloc(num, PS_TYPE_F32); + psVector *refBG = psVectorAlloc(num, PS_TYPE_F32); + + psVector *inFlux = psVectorAlloc(num, PS_TYPE_F32); + psVector *inBG = psVectorAlloc(num, PS_TYPE_F32); + psMinimization *min = psMinimizationAlloc(MAXITER, TOLERANCE); + for (int i = 0; i < num; i++) { + psPlane *coords = starList->data[i]; + + // Iterate over postage stamp + int xCen = (int)(coords->x + 0.5); + int yCen = (int)(coords->y + 0.5); + int xMin = MAX(0, xCen - size); + int xMax = MIN(refImage->numCols - 1, xCen + size); + int yMin = MAX(0, yCen - size); + int yMax = MIN(refImage->numRows - 1, yCen + size); + + // Tweak position if necessary + if (xCen >= refImage->numCols) { + xCen = refImage->numCols - 1; + } + if (yCen >= refImage->numRows) { + yCen = refImage->numRows - 1; + } + if (xCen < 0) { + xCen = 0; + } + if (yCen < 0) { + yCen = 0; + } + + psTrace("phot", 3, "Postage stamp: %d,%d-->%d,%d\n", xMin, yMin, xMax, yMax); + + int numPix = (xMax - xMin + 1) * (yMax - yMin + 1); + psArray *pixels = psArrayAlloc(numPix); + psVector *fluxes = psVectorAlloc(numPix, PS_TYPE_F32); + psVector *errors = psVectorAlloc(numPix, PS_TYPE_F32); + + int index = 0; + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + psVector *pixel = psVectorAlloc(2, PS_TYPE_F32); + pixel->data.F32[0] = x; + pixel->data.F32[1] = y; + pixels->data[index] = pixel; + fluxes->data.F32[index] = refImage->data.F32[y][x]; + errors->data.F32[index] = sqrt(refImage->data.F32[y][x]); + index++; + } + } + + psVector *params = psVectorAlloc(7, PS_TYPE_F32); + params->data.F32[0] = refMedian;// Background + params->data.F32[1] = refImage->data.F32[yCen][xCen]; // Normalisation + params->data.F32[2] = coords->x; // x position + params->data.F32[3] = coords->y; // y position + params->data.F32[4] = 1.0; // width in x + params->data.F32[5] = 1.0; // width in y + params->data.F32[6] = M_PI/4; // width in xy + + bool fitOK = psMinimizeLMChi2(min, NULL, params, NULL, pixels, fluxes, NULL, gaussian2d_deriv); + //bool fitOK = psMinimizeChi2Powell(min, params, NULL, pixels, fluxes, NULL, gaussian2d_noderiv); + + printf("Final result: %d %d %f --> %f %f %f %f %f %f %f\n", fitOK, min->iter, min->value, + params->data.F32[0], params->data.F32[1], params->data.F32[2], params->data.F32[3], + params->data.F32[4], params->data.F32[5], params->data.F32[6]); + + refFlux->data.F32[i] = params->data.F32[1] * 2.0 * M_PI * params->data.F32[4] * params->data.F32[5]; + refBG->data.F32[i] = params->data.F32[0]; + + printf("%f %f --> %f %f %f %f\n", coords->x, coords->y, params->data.F32[2], params->data.F32[3], + refFlux->data.F32[i], refBG->data.F32[i]); + + } + + for (int i = 0; i < num; i++) { + psPlane *coords = starList->data[i]; + printf("%f %f --> %f %f\n", coords->x, coords->y, -2.5*log10f(refFlux->data.F32[i]) + ZP, + refBG->data.F32[i]); + } + +} Index: /trunk/stac/src/shift.c =================================================================== --- /trunk/stac/src/shift.c (revision 3613) +++ /trunk/stac/src/shift.c (revision 3613) @@ -0,0 +1,100 @@ +#include +#include "pslib.h" +#include "stac.h" +#include + +double getTime(void) +/* Gets the current time. Got this from Nick Kaiser's fetchpix.c */ +{ + struct timeval tv; + gettimeofday(&tv, NULL); + return(tv.tv_sec + 1.e-6 * tv.tv_usec); +} + + +int main(int argc, char **argv) +{ + + double startTime = getTime(); + +#if 0 + psMemAllocateCallbackSet(stacMemPrint); + psMemAllocateCallbackSetID(185); + psMemFreeCallbackSet(stacMemPrint); + psMemFreeCallbackSetID(185); +#endif + + // Set trace levels + psTraceSetLevel(".",0); + psTraceSetLevel("stac.checkMemory",10); + psTraceSetLevel("stac.config",10); + psTraceSetLevel("stac.read",10); + psTraceSetLevel("stac.invertMaps",10); + psTraceSetLevel("stac.errors",10); + psTraceSetLevel("stac.transform",10); + psTraceSetLevel("stac.combine",10); + psTraceSetLevel("stac.rejection",10); + psTraceSetLevel("stac.time",10); + psTraceSetLevel("stac.area",10); + psTraceSetLevel("stac.size",10); + + // Set logging level + psLogSetLevel(9); + + // Parse command line + stacConfig *config = stacParseConfig(argc, argv); + + // Read input files + psArray *inputs = stacReadImages(config); + + // Read maps + psArray *maps = stacReadMaps(config); + + psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime); + + // Get size + stacSize(config, inputs, maps); + + // Invert maps + psArray *inverseMaps = stacInvertMaps(maps, config); + + // Generate errors + psArray *errors = stacErrorImages(inputs, config); + + // Transform inputs and errors + psArray *transformed = NULL; + psArray *transformedErrors = NULL; + (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config); + + psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime); + + // Combine with rejection + psArray *rejected = NULL; + psImage *combined = NULL; + (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config); + + psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime); + + psFits *outFile = psFitsAlloc(config->output); + if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); + } + psTrace("stac", 1, "Combined image written to %s\n", config->output); + psFree(outFile); + + // Free everything I've used + stacConfigFree(config); + psFree(inputs); + psFree(maps); + psFree(inverseMaps); + psFree(errors); + psFree(transformedErrors); + psFree(transformed); + psFree(rejected); + psFree(combined); + + psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime); + + // Check memory for leaks, corruption + stacCheckMemory(); +} Index: /trunk/stac/src/stacScales.c =================================================================== --- /trunk/stac/src/stacScales.c (revision 3613) +++ /trunk/stac/src/stacScales.c (revision 3613) @@ -0,0 +1,263 @@ +#include +#include +#include "pslib.h" +#include "stac.h" + +#define SQUARE(x) ((x)*(x)) + +#define SAMPLE 10 // Subsample rate for images + +float stacBackground(const psImage *image, // Image for which to get the background + int sample // Sample in increments of this value + ) +{ + assert(image->type.type == PS_TYPE_F32); + +// Will use robust median instead of sampling --- it's supposed to be fast. +#if 1 + int size = image->numCols * image->numRows; // Number of pixels in image + int numSamples = size / sample; // Number of samples in image + psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample + + int offset = 0; // Offset from start of the row + int index = 0; // Sample number + for (int row = 0; row < image->numRows; row++) { + // I'd cast this as a "for", but this makes it a bit easier to understand. + int col = offset; + while (col < image->numCols) { + values->data.F32[index] = image->data.F32[row][col]; + col += sample; + index++; + } + offset = col - image->numCols; + } + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + stats = psVectorStats(stats, values, NULL, NULL, 0); + float median = stats->sampleMedian; + psFree(stats); + psFree(values); +#else + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Using a clipped mean because median is SLOW + stats->clipSigma = 3.0; + stats->clipIter = 5; + stats = psImageStats(stats, (psImage*)image, NULL, 0); + float median = stats->robustMedian; + psFree(stats); +#endif + + return median; +} + + +bool stacScales(psVector **scalesPtr, // Scales to return + psVector **offsetsPtr, // Offsets to return + const psArray *images, // Images on which to measure the scales and offsets + const stacConfig *config // Configuration + ) +{ + assert(scalesPtr); + assert(offsetsPtr); + assert(images); + assert(config); + for (int i = 0; i < images->n; i++) { + psImage *image = images->data[i]; + assert(image->type.type == PS_TYPE_F32); + } + + psVector *scales = NULL; // Relative scales between images + psVector *offsets = NULL; // Offsets between images (ADU) + if (*scalesPtr) { + scales = *scalesPtr; + assert(scales); + assert(scales->n == images->n); + assert(scales->type.type == PS_TYPE_F32); + } else { + scales = psVectorAlloc(images->n, PS_TYPE_F32); + *scalesPtr = scales; + } + if (*offsetsPtr) { + offsets = *offsetsPtr; + assert(offsets); + assert(offsets->n == images->n); + assert(offsets->type.type == PS_TYPE_F32); + } else { + offsets = psVectorAlloc(images->n, PS_TYPE_F32); + *offsetsPtr = offsets; + } + + // Offsets are easy + psTrace("stac.scales", 3, "Getting background levels....\n"); + double startTime = getTime(); + for (int i = 0; i < images->n; i++) { + offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE); + psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]); + double time = getTime(); + psTrace("stac.scales", 10, "Took %f sec\n", time - startTime); + startTime = time; + } + + // Now the scales + if (config->starFile == NULL || config->starMapFile == NULL) { + psLogMsg("stac.scales", PS_LOG_INFO, + "No coordinates available to set scales --- assuming all are identical.\n"); + for (int i = 0; i < images->n; i++) { + scales->data.F32[i] = 1.0; + psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]); + } + } else { + // Read star coordinates and map + psArray *starCoords = stacReadCoords(config->starFile); // Array of star coordinates + psPlaneTransform *starMap = stacReadMap(config->starMapFile); // Transformation for star coordinates + + // Transform the stellar positions to match the transformed reference frame + psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions + // Fix up difference between map and output frame + starMap->x->coeff[0][0] -= config->xMapDiff; + starMap->y->coeff[0][0] -= config->yMapDiff; + for (int i = 0; i < starCoords->n; i++) { + starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]); + } + + psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors + psArray *masks = psArrayAlloc(images->n); // Array of masks for stars + + // Set scales relative to the first image + scales->data.F32[0] = 1.0; + + // Iterate over images + for (int i = 0; i < scales->n; i++) { + psImage *image = images->data[i]; // The image we're working with + + float background = offsets->data.F32[i]; // Background in image + + // Do photometry on transformed image i + psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i); + psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars + psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry + for (int j = 0; j < starCoordsTransformed->n; j++) { + psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star + + if (coords->x < config->aper || coords->y < config->aper || + coords->x + config->aper > image->numCols - 1 || + coords->y + config->aper > image->numRows) { + mask->data.U8[j] = 1; + } else { + // Sum flux within the aperture + float sum = 0.0; + int numPix = 0; + float aper2 = SQUARE(config->aper); + for (int y = (int)floorf(coords->y - config->aper); + y <= (int)ceilf(coords->y + config->aper); y++) { + for (int x = (int)floorf(coords->x - config->aper); + x <= (int)ceilf(coords->x + config->aper); x++) { + if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <= + aper2) { + sum += image->data.F32[y][x]; + numPix++; + } + } + } + psTrace("stac.scales", 9, "Star at %f,%f --> %f\n", coords->x, coords->y, sum); + // Subtract background, renormalise to account for circular aperture + if (numPix > 0 && sum > 0) { + sum -= offsets->data.F32[i] * (float)numPix; + photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix; + mask->data.U8[j] = 0; + } else { + mask->data.U8[j] = 1; + } + } + } + stars->data[i] = photometry; + masks->data[i] = mask; + } + + // Get the scales + psVector *ref = stars->data[0]; // The reference photometry + psVector *refMask = masks->data[0]; // The reference mask + psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics + stats->clipSigma = 2.5; + stats->clipIter = 3; + for (int i = 1; i < scales->n; i++) { + psVector *compare = stars->data[i]; // The comparison photometry + psVector *compareMask = masks->data[i]; // The comparison mask + + compare = (psVector*)psBinaryOp(compare, compare, "/", ref); + compareMask = (psVector*)psBinaryOp(compareMask, compareMask, "+", refMask); + + stats = psVectorStats(stats, compare, NULL, compareMask, 3); // Use maskVal of 3 to catch 1 and 2 + + scales->data.F32[i] = stats->clippedMean; + psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]); + } + psFree(stats); + + psFree(stars); + psFree(masks); + } + + // Change the saturation and bad values + psVector *saturated = config->saturated; // Saturation limits + psVector *bad = config->bad; // Bad limits + for (int i = 0; i < saturated->n; i++) { + saturated->data.F32[i] = (saturated->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i]; + bad->data.F32[i] = (bad->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i]; + } + + return true; +} + + +bool stacRescale(psArray *images, // Images to rescale + psArray *errImages, // Variance images to rescale + const psImage *mask, // Mask for pixels to scale + const psVector *scales,// Scales for images + const psVector *offsets // Offsets for images + ) +{ + assert(images); + assert(scales); + assert(offsets); + assert(images->n == scales->n); + assert(images->n == offsets->n); + assert(!errImages || errImages->n == images->n); + assert(scales->type.type == PS_TYPE_F32); + assert(offsets->type.type == PS_TYPE_F32); + for (int i = 0; i < images->n; i++) { + psImage *image = images->data[i]; // Image of interest + assert(image->type.type == PS_TYPE_F32); + assert(scales->data.F32[i] != 0); + if (mask) { + assert(mask->type.type == PS_TYPE_U8); + assert(image->numCols == mask->numCols && image->numRows == mask->numRows); + } + if (errImages) { + psImage *errImage = errImages->data[i]; + assert(errImage->type.type == PS_TYPE_F32); + assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows); + } + } + + for (int i = 0; i < images->n; i++) { + psImage *image = images->data[i]; // Image to rescale + psImage *errImage = NULL; // Variance image to rescale + if (errImages) { + errImage = errImages->data[i]; + } + float scale = scales->data.F32[i]; // Scale to use + float offset = offsets->data.F32[i]; // Offset to use + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + if (!mask || mask->data.F32[y][x]) { + image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale; + if (errImage) { + errImage->data.F32[y][x] /= SQUARE(scale); + } + } + } + } + } + + return true; +}