Index: /trunk/stac/src/Makefile =================================================================== --- /trunk/stac/src/Makefile (revision 2763) +++ /trunk/stac/src/Makefile (revision 2764) @@ -1,10 +1,10 @@ SHELL = /bin/sh CC = gcc -CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include # -DCRFLUX -DTESTING +CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX -DTESTING PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm LDFLAGS += $(PSLIB) OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \ - stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o + stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o TARGET = stac Index: /trunk/stac/src/stac.c =================================================================== --- /trunk/stac/src/stac.c (revision 2763) +++ /trunk/stac/src/stac.c (revision 2764) @@ -72,6 +72,20 @@ #endif + // Get regions of interest in the source frame + psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest + for (int i = 0; i < inputs->n; i++) { + fprintf(stderr, "Finding area of interest for image %d\n", i); + regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i], + ((psImage*)(inputs->data[i]))->numCols, + ((psImage*)(inputs->data[i]))->numRows); +#ifdef TESTING + char regionFile[MAXCHAR]; // Filename of region image + sprintf(regionFile,"%s.region",config->inputs->data[i]); + psImageWriteSection(regions->data[i], 0, 0, 0, NULL, 0, regionFile); +#endif + } + // Transform rejected pixels to source frame - psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, inverseMaps, rejected, config); + psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config); // Redo transformation with the mask @@ -93,4 +107,5 @@ // Free everything I've used stacConfigFree(config); + psFree(regions); psFree(inputs); psFree(maps); Index: /trunk/stac/src/stac.h =================================================================== --- /trunk/stac/src/stac.h (revision 2763) +++ /trunk/stac/src/stac.h (revision 2764) @@ -137,12 +137,29 @@ // Transform the rejection masks back to the source frame psArray *stacRejection(psArray *inputs, // Input images - psArray *transformed, // Transformed images - psImage *combined, // Combined image + psArray *rejected, // Rejected images + psArray *regions, // Regions of interest psArray *maps, // Maps from input to transformed image psArray *inverseMaps, // Maps from transformed to input image - psArray *rejected, // Rejected images stacConfig *config // Configuration ); +/************************************************************************************************************/ + +// stacAreaOfInterest.c + +// Returns the change in x' and y' for a change in x and y of 1 +bool stacPlaneTransformDeriv(psPlane *out, // Output derivative, assumed already allocated + psPlaneTransform *transform, // The transform for which to obtain the derivative + psPlane *in // The position at which to obtain the derivative + ); + +// Return a mask image designating the region of interest in the source frame +// Basically, uses derivatives of the transformation to work out what pixels in the source might be of interest +psImage *stacAreaOfInterest(psImage *mask, // Mask that is to be used to determine the area of interest + psPlaneTransform *map, // Map from the mask to the area of interest + int nxOut, int nyOut // Size of the area of interest + ); + +/************************************************************************************************************/ #endif Index: /trunk/stac/src/stacAreaOfInterest.c =================================================================== --- /trunk/stac/src/stacAreaOfInterest.c (revision 2764) +++ /trunk/stac/src/stacAreaOfInterest.c (revision 2764) @@ -0,0 +1,125 @@ +#include +#include "pslib.h" +#include "stac.h" + +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#define MIN(a,b) ((a) < (b) ? (a) : (b)) + +// Returns the change in x' and y' for a change in x and y of 1 +bool stacPlaneTransformDeriv(psPlane *out, // Output derivative, assumed already allocated + psPlaneTransform *transform, // The transform for which to obtain the derivative + psPlane *in // The position at which to obtain the derivative + ) +{ + psDPolynomial2D *xPoly = transform->x; + psDPolynomial2D *yPoly = transform->y; + + if (xPoly->type != PS_POLYNOMIAL_ORD || yPoly->type != PS_POLYNOMIAL_ORD) { + psError("stac.area.deriv", "Transformation polynomials aren't ordinary polynomials.\n"); + return false; + } + + // Initialise derivatives + out->x = 0.0; + out->y = 0.0; + + // Do x component + double xSum = 1.0 / in->x; // 1/x in order to calculate x^(i-1) + for (int xOrder = 0; xOrder < xPoly->nX; xOrder++) { + double ySum = xSum / in->y; // 1/y in order to calculate y^(j-1) + for (int yOrder = 0; yOrder < xPoly->nY; yOrder++) { + if (xPoly->mask[xOrder][yOrder] == 0) { + out->x += xPoly->coeff[xOrder][yOrder] * (double)xOrder * ySum * in->y; // a_ij i x^(i-1) y^j + out->x += xPoly->coeff[xOrder][yOrder] * (double)yOrder * ySum * in->x; // a_ij j x^i y^(j-1) + } + ySum *= in->y; + } + xSum *= in->x; + } + + // Do y component + xSum = 1.0 / in->x; // 1/x in order to calculate x^(i-1) + for (int xOrder = 0; xOrder < yPoly->nX; xOrder++) { + double ySum = xSum / in->y; // 1/y in order to calculate y^(j-1) + for (int yOrder = 0; yOrder < yPoly->nY; yOrder++) { + if (yPoly->mask[xOrder][yOrder] == 0) { + out->y += yPoly->coeff[xOrder][yOrder] * (float)xOrder * ySum * in->y; // a_ij i x^(i-1) y^j + out->y += yPoly->coeff[xOrder][yOrder] * (float)yOrder * ySum * in->x; // a_ij j x^i y^(j-1) + } + ySum *= in->y; + } + xSum *= in->x; + } + + return true; +} + + + +// Calculate +psImage *stacAreaOfInterest(psImage *mask, // Mask that is to be used to determine the area of interest + psPlaneTransform *map, // Map from the mask to the area of interest + int nxOut, int nyOut // Size of the area of interest + ) +{ + // Input checks + if (mask == NULL) { + psError("stac.area", "Input mask is NULL for some weird reason.\n"); + return NULL; + } + if (map == NULL) { + psError("stac.area", "Input map is NULL for some weird reason.\n"); + return NULL; + } + if (mask->type.type != PS_TYPE_U8) { + psError("stac.area", "Input mask must be of type U8.\n"); + return NULL; + } + + // Number of pixels in x and y + int nxIn = mask->numCols; + int nyIn = mask->numRows; + + // Create and zero output image + psImage *output = psImageAlloc(nxOut, nyOut, PS_TYPE_U8); + for (int y = 0; y < nyOut; y++) { + for (int x = 0; x < nxOut; x++) { + output->data.U8[y][x] = 0; + } + } + + // Coordinates for the transformations + psPlane *inFrame = psAlloc(sizeof(psPlane)); + psPlane *outFrame = psAlloc(sizeof(psPlane)); + psPlane *deriv = psAlloc(sizeof(psPlane)); + + // Find pixels on the mask that are of interest + for (int y = 0; y < nyIn; y++) { + for (int x = 0; x < nxIn; x++) { + if (mask->data.U8[y][x]) { + inFrame->x = (double)x + 0.5; + inFrame->y = (double)y + 0.5; + (void)psPlaneTransformApply(outFrame, map, inFrame); + (void)stacPlaneTransformDeriv(deriv, map, inFrame); + + // Calculate range + int iMin = MAX((int)(outFrame->x - 0.5*deriv->x), 0); + int iMax = MIN((int)(outFrame->x + 0.5*deriv->x + 0.5), nxOut - 1); + int jMin = MAX((int)(outFrame->y - 0.5*deriv->y), 0); + int jMax = MIN((int)(outFrame->y + 0.5*deriv->y + 0.5), nyOut - 1); + + for (int j = jMin; j <= jMax; j++) { + for (int i = iMin; i <= iMax; i++) { + output->data.U8[j][i] = 1; + } + } + } // Pixel of interest + } + } // Iterating over input mask + + psFree(inFrame); + psFree(outFrame); + psFree(deriv); + + return output; +} Index: /trunk/stac/src/stacCombine.c =================================================================== --- /trunk/stac/src/stacCombine.c (revision 2763) +++ /trunk/stac/src/stacCombine.c (revision 2764) @@ -106,9 +106,9 @@ } for (int i = 0; i < nImages; i++) { - (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Create rejection mask + (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8); // Create rejection mask // Zero out the mask for (int r = 0; r < numRows; r++) { for (int c = 0; c < numCols; c++) { - ((psImage*)((*rejected)->data[i]))->data.F32[r][c] = 0.0; + ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0; } } @@ -154,7 +154,8 @@ // Rejection iterations - for (int rejNum = 0; rejNum < nReject; rejNum++) { - float max = 0.0; - int maxIndex = -1; + bool keepGoing = true; // Keep going with rejection? + for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) { + float max = 0.0; // Maximum deviation + int maxIndex = -1; // Index of the maximum deviation for (int i = 0; i < nImages; i++) { if (!mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) { @@ -164,13 +165,15 @@ } // Reject the pixel with the maximum deviation - if (max > reject) { + if (max > reject) { mask->data.U8[maxIndex] = 1; - ((psImage*)((*rejected)->data[maxIndex]))->data.F32[y][x] += 1.0; + ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1.0; // Re-do combination following rejection average = stacCombineMean(pixels, deltas, mask); - } - } // Rejection iterations - - combined->data.F32[y][x] = stacCombineMean(pixels, deltas, mask); + } else { + keepGoing = false; + } + } + + combined->data.F32[y][x] = average; } } // Iterating over output pixels Index: /trunk/stac/src/stacRead.c =================================================================== --- /trunk/stac/src/stacRead.c (revision 2763) +++ /trunk/stac/src/stacRead.c (revision 2764) @@ -81,6 +81,6 @@ // Mask out the cross-terms - map->x->mask[1][1] = 0; - map->y->coeff[1][1] = 0; + map->x->mask[1][1] = 1; + map->y->mask[1][1] = 1; // Set the cross-terms to zero, just in case... map->x->coeff[1][1] = 0.0; Index: /trunk/stac/src/stacRejection.c =================================================================== --- /trunk/stac/src/stacRejection.c (revision 2763) +++ /trunk/stac/src/stacRejection.c (revision 2764) @@ -42,41 +42,9 @@ } -#if 0 -float stacGradient(psImage *image, // Input for which to measure the gradient - int x, int y // Coordinates at which to measure the gradient - ) -{ - float sum = 0.0; // The sum of surrounding pixels - float maxDiff = 0.0; - int num = 0; - // Get limits - int xMin = MAX(x - 1, 0); - int xMax = MIN(x + 1, image->numCols - 1); - int yMin = MAX(y - 1, 0); - int yMax = MIN(y + 1, image->numRows - 1); - for (int j = yMin; j <= yMax; j++) { - for (int i = xMin; i <= xMax; i++) { - if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) { - maxDiff = image->data.F32[j][i] - image->data.F32[y][x]; - } - sum += image->data.F32[j][i]; - num++; - } - } - sum -= image->data.F32[y][x]; - sum /= (float)(num-1); - sum -= image->data.F32[y][x]; -// return sum / image->data.F32[y][x]; - return maxDiff / image->data.F32[y][x]; -} -#endif - - psArray *stacRejection(psArray *inputs, // Input images - psArray *transformed, // Transformed images - psImage *combined, // Combined image + psArray *rejected, // Rejected images + psArray *regions, // Regions of interest psArray *maps, // Maps from input to transformed image psArray *inverseMaps, // Maps from transformed to input image - psArray *rejected, // Rejected images stacConfig *config // Configuration ) @@ -85,31 +53,34 @@ // Check inputs - if (inputs->n != transformed->n) { - psError("stac.rejection", "Number of input (%d) and transformed (%d) images does not match.\n", - inputs->n, transformed->n); + if (inputs->n != rejected->n) { + psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n", + inputs->n, rejected->n); + return NULL; + } + if (inputs->n != regions->n) { + psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n", + inputs->n, regions->n); + return NULL; + } + if (inputs->n != maps->n) { + psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does " + "not match.\n", inputs->n, maps->n); + return NULL; + } + if (inputs->n != inverseMaps->n) { + psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does " + "not match.\n", inputs->n, inverseMaps->n); return NULL; } for (int i = 0; i < nImages; i++) { - if ((((psImage*)(transformed->data[i]))->numRows != combined->numRows) || - (((psImage*)(transformed->data[i]))->numCols != combined->numCols)) { + if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) || + (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) { psError("stac.rejection", - "Sizes of transformed (%dx%d) and combined (%dx%d) images do not match.\n", - ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows, - combined->numCols, combined->numRows); + "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n", + ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows, + ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i); return NULL; } - if ((((psImage*)(transformed->data[i]))->numRows != ((psImage*)(rejected->data[i]))->numRows) || - (((psImage*)(transformed->data[i]))->numCols != ((psImage*)(rejected->data[i]))->numCols)) { - psError("stac.rejection", - "Sizes of transformed (%dx%d) and rejected (%dx%d) images do not match.\n", - ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows, - ((psImage*)(rejected->data[i]))->numCols, ((psImage*)(rejected->data[i]))->numRows); - return NULL; - } - } - - // Size of output images - int nxOutput = combined->numCols; - int nyOutput = combined->numRows; + } // Stuff for the transformations @@ -132,4 +103,5 @@ psImage *reject = rejected->data[i]; // Pull out the mask in the output frame psPlaneTransform *map = maps->data[i]; // The map from input to output + psImage *region = regions->data[i]; // The region of interest for this image psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i); @@ -153,45 +125,40 @@ for (int y = 0; y < nyInput; y++) { for (int x = 0; x < nxInput; x++) { - inCoords->x = (double)x; - inCoords->y = (double)y; - (void)psPlaneTransformApply(outCoords, map, inCoords); - float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y, - NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + + // Only transform pixels of interest + if (region->data.U8[y][x]) { + + inCoords->x = (double)x; + inCoords->y = (double)y; + (void)psPlaneTransformApply(outCoords, map, inCoords); + float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y, + NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); #ifdef TESTING - rejmap->data.F32[y][x] = maskVal; -#endif - -#ifdef TESTING - // Check gradient - float meanGrads = 0.0; - int numGrads = 0; - for (int j = 0; j < nImages; j++) { - if (i != j) { - // Get coordinates for that image - (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords); - int xPix = (int)(inCoords->x + 0.5); - int yPix = (int)(inCoords->y + 0.5); - if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) && - (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) { - // Calculate the gradient - meanGrads += stacGradient(inputs->data[j], xPix, yPix); - numGrads++; + rejmap->data.F32[y][x] = maskVal; + + // Calculate gradient + float meanGrads = 0.0; + int numGrads = 0; + for (int j = 0; j < nImages; j++) { + if (i != j) { + // Get coordinates for that image + (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords); + int xPix = (int)(inCoords->x + 0.5); + int yPix = (int)(inCoords->y + 0.5); + if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) && + (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) { + // Calculate the gradient + meanGrads += stacGradient(inputs->data[j], xPix, yPix); + numGrads++; + } } } - } - if (numGrads > 0) { - meanGrads /= (float)numGrads; - } else { - meanGrads = 0; - } -#else - float meanGrads = 10000.0; -#endif - - -#ifdef TESTING - grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads; -#endif - if (maskVal > config->frac) { + if (numGrads > 0) { + meanGrads /= (float)numGrads; + } else { + meanGrads = 0; + } + grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads; +#endif if ((maskVal > config->frac) && @@ -200,11 +167,13 @@ nBad++; #ifdef CRFLUX - fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x], - maskVal, stacGradient(inputs->data[i], x, y)); -#endif + fprintf(crs, "%d %d --> %f %f %f\n", x, y, + ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal, + stacGradient(inputs->data[i], x, y)); +#endif + } else { + mask->data.U8[y][x] = 0; } - } else { - mask->data.U8[y][x] = 0; - } + + } // Only touching pixels of interest }