Changeset 5742 for trunk/pois/src/poisCalculateDeviations.c
- Timestamp:
- Dec 7, 2005, 2:35:14 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/pois/src/poisCalculateDeviations.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/pois/src/poisCalculateDeviations.c
r5717 r5742 6 6 #define MAXCHAR 80 7 7 8 psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL9 psArray *stamps, // Array of stamps10 psImage *refImage, // Reference image11 psImage *inImage, // Input image12 psImage *mask, // Mask image13 psArray *kernelParams, // Array of kernel parameters14 psVector *solution, // Solution vector15 poisConfig *config // Configuration8 psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL 9 psArray *stamps, // Array of stamps 10 psImage *refImage, // Reference image 11 psImage *inImage, // Input image 12 psImage *mask, // Mask image 13 psArray *kernelParams, // Array of kernel parameters 14 psVector *solution, // Solution vector 15 poisConfig *config // Configuration 16 16 ) 17 17 { … … 32 32 33 33 if (!deviations) { 34 deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);34 deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); 35 35 } 36 36 37 int footprint = config->footprint; // Size of stamp footprint37 int footprint = config->footprint; // Size of stamp footprint 38 38 int xSize = footprint + config->xKernel; // (Half) size of subimage in x 39 39 int ySize = footprint + config->yKernel; // (Half) size of subimage in y 40 psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics40 psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics 41 41 42 42 psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp 43 43 for (int s = 0; s < stamps->n; s++) { 44 poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest45 int x = stamp->x;// Stamp x coord46 int y = stamp->y;// Stamp y coord47 if (stamp->status == POIS_STAMP_USED) {48 psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};49 psImage *refStamp = psImageSubset(refImage, stampRegion);50 psImage *inStamp = psImageSubset(inImage, stampRegion);51 psImage *maskStamp = psImageSubset(mask, stampRegion);52 psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);53 // Calculate chi^254 (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);55 (void)psBinaryOp(subStamp, subStamp, "*", subStamp);56 (void)psBinaryOp(subStamp, subStamp, "/", inStamp);57 psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,58 config->yKernel + 2 * footprint };59 psImage *subStampTrim = psImageSubset(subStamp, stampTrim);60 psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);61 // Copy image to workaround bug 30562 psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);63 psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);64 65 (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);66 67 psFree(tempImage);68 psFree(tempMask);69 70 deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);71 psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,72 deviations->data.F32[s]);73 44 poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest 45 int x = stamp->x; // Stamp x coord 46 int y = stamp->y; // Stamp y coord 47 if (stamp->status == POIS_STAMP_USED) { 48 psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize}; 49 psImage *refStamp = psImageSubset(refImage, stampRegion); 50 psImage *inStamp = psImageSubset(inImage, stampRegion); 51 psImage *maskStamp = psImageSubset(mask, stampRegion); 52 psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config); 53 // Calculate chi^2 54 (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp); 55 (void)psBinaryOp(subStamp, subStamp, "*", subStamp); 56 (void)psBinaryOp(subStamp, subStamp, "/", inStamp); 57 psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel, 58 config->yKernel + 2 * footprint }; 59 psImage *subStampTrim = psImageSubset(subStamp, stampTrim); 60 psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim); 61 // Copy image to workaround bug 305 62 psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32); 63 psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8); 64 65 (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD); 66 67 psFree(tempImage); 68 psFree(tempMask); 69 70 deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0); 71 psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y, 72 deviations->data.F32[s]); 73 74 74 #ifdef TESTING 75 char stampName[MAXCHAR];// File name for stamp76 snprintf(stampName, MAXCHAR, "stamp%d.fits", s);77 psFits *stampFile = psFitsAlloc(stampName);78 if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {79 psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);80 }81 psFree(stampFile);75 char stampName[MAXCHAR]; // File name for stamp 76 snprintf(stampName, MAXCHAR, "stamp%d.fits", s); 77 psFits *stampFile = psFitsOpen(stampName, "w"); 78 if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) { 79 psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName); 80 } 81 psFitsClose(stampFile); 82 82 #endif 83 83 84 84 #if 0 85 psFree(convRefStamp);86 psFree(maskStampTrim);87 psFree(subStampTrim);88 psFree(maskStamp);89 psFree(refStamp);90 psFree(inStamp);85 psFree(convRefStamp); 86 psFree(maskStampTrim); 87 psFree(subStampTrim); 88 psFree(maskStamp); 89 psFree(refStamp); 90 psFree(inStamp); 91 91 #endif 92 }92 } 93 93 } 94 94
Note:
See TracChangeset
for help on using the changeset viewer.
