Changeset 21211
- Timestamp:
- Jan 28, 2009, 2:33:51 PM (17 years ago)
- Location:
- branches/pap_branch_20090128/psModules/src
- Files:
-
- 41 edited
-
camera/pmFPA.c (modified) (9 diffs)
-
camera/pmFPA.h (modified) (6 diffs)
-
camera/pmFPAConstruct.c (modified) (3 diffs)
-
camera/pmFPACopy.c (modified) (3 diffs)
-
camera/pmFPAExtent.c (modified) (4 diffs)
-
camera/pmFPAMaskWeight.c (modified) (37 diffs)
-
camera/pmFPAMaskWeight.h (modified) (7 diffs)
-
camera/pmFPAMosaic.c (modified) (24 diffs)
-
camera/pmFPARead.c (modified) (15 diffs)
-
camera/pmFPARead.h (modified) (8 diffs)
-
camera/pmFPAWrite.c (modified) (7 diffs)
-
camera/pmFPAWrite.h (modified) (9 diffs)
-
camera/pmFPAfile.c (modified) (2 diffs)
-
camera/pmFPAfile.h (modified) (2 diffs)
-
camera/pmFPAfileFitsIO.c (modified) (4 diffs)
-
camera/pmFPAfileFitsIO.h (modified) (3 diffs)
-
camera/pmFPAfileIO.c (modified) (10 diffs)
-
camera/pmHDU.c (modified) (5 diffs)
-
camera/pmHDU.h (modified) (5 diffs)
-
camera/pmHDUGenerate.c (modified) (11 diffs)
-
camera/pmHDUUtils.c (modified) (2 diffs)
-
camera/pmReadoutStack.c (modified) (7 diffs)
-
detrend/pmShutterCorrection.c (modified) (8 diffs)
-
imcombine/pmPSFEnvelope.c (modified) (5 diffs)
-
imcombine/pmReadoutCombine.c (modified) (14 diffs)
-
imcombine/pmReadoutCombine.h (modified) (4 diffs)
-
imcombine/pmStack.c (modified) (6 diffs)
-
imcombine/pmSubtraction.c (modified) (17 diffs)
-
imcombine/pmSubtractionEquation.c (modified) (25 diffs)
-
imcombine/pmSubtractionMatch.c (modified) (9 diffs)
-
imcombine/pmSubtractionParams.c (modified) (10 diffs)
-
imcombine/pmSubtractionStamps.c (modified) (7 diffs)
-
imcombine/pmSubtractionStamps.h (modified) (2 diffs)
-
imcombine/pmSubtractionThreads.c (modified) (2 diffs)
-
objects/pmSource.c (modified) (14 diffs)
-
objects/pmSource.h (modified) (3 diffs)
-
objects/pmSourceFitModel.c (modified) (4 diffs)
-
objects/pmSourceFitSet.c (modified) (15 diffs)
-
objects/pmSourceMoments.c (modified) (12 diffs)
-
objects/pmSourcePhotometry.c (modified) (10 diffs)
-
objects/pmSourceSky.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/pap_branch_20090128/psModules/src/camera/pmFPA.c
r19013 r21211 32 32 psFree(readout->image); 33 33 psFree(readout->mask); 34 psFree(readout->weight); 34 psFree(readout->variance); 35 psFree(readout->covariance); 35 36 psFree(readout->analysis); 36 37 psFree(readout->bias); … … 163 164 psFree(readout->image); 164 165 psFree(readout->mask); 165 psFree(readout->weight); 166 psFree(readout->variance); 167 psFree(readout->covariance); 166 168 psFree(readout->bias); 167 169 … … 169 171 170 172 readout->image = NULL; 171 readout->weight = NULL; 173 readout->variance = NULL; 174 readout->covariance = NULL; 172 175 readout->mask = NULL; 173 176 … … 179 182 readout->thisImageScan = 0; 180 183 readout->thisMaskScan = 0; 181 readout->this WeightScan = 0;184 readout->thisVarianceScan = 0; 182 185 183 186 readout->lastImageScan = 0; 184 187 readout->lastMaskScan = 0; 185 readout->last WeightScan = 0;188 readout->lastVarianceScan = 0; 186 189 } 187 190 … … 197 200 if (cell->hdu) { 198 201 psFree(cell->hdu->images); 199 psFree(cell->hdu-> weights);202 psFree(cell->hdu->variances); 200 203 psFree(cell->hdu->masks); 201 204 // psFree(cell->hdu->header); 202 205 203 206 cell->hdu->images = NULL; 204 cell->hdu-> weights = NULL;207 cell->hdu->variances = NULL; 205 208 cell->hdu->masks = NULL; 206 209 // cell->hdu->header = NULL; … … 219 222 if (chip->hdu) { 220 223 psFree(chip->hdu->images); 221 psFree(chip->hdu-> weights);224 psFree(chip->hdu->variances); 222 225 psFree(chip->hdu->masks); 223 226 // psFree(chip->hdu->header); 224 227 225 228 chip->hdu->images = NULL; 226 chip->hdu-> weights = NULL;229 chip->hdu->variances = NULL; 227 230 chip->hdu->masks = NULL; 228 231 // chip->hdu->header = NULL; … … 241 244 if (fpa->hdu) { 242 245 psFree(fpa->hdu->images); 243 psFree(fpa->hdu-> weights);246 psFree(fpa->hdu->variances); 244 247 psFree(fpa->hdu->masks); 245 248 // psFree(fpa->hdu->header); 246 249 247 250 fpa->hdu->images = NULL; 248 fpa->hdu-> weights = NULL;251 fpa->hdu->variances = NULL; 249 252 fpa->hdu->masks = NULL; 250 253 // fpa->hdu->header = NULL; … … 259 262 tmpReadout->image = NULL; 260 263 tmpReadout->mask = NULL; 261 tmpReadout->weight = NULL; 264 tmpReadout->variance = NULL; 265 tmpReadout->covariance = NULL; 262 266 tmpReadout->bias = psListAlloc(NULL); 263 267 tmpReadout->analysis = psMetadataAlloc(); … … 276 280 tmpReadout->thisImageScan = 0; 277 281 tmpReadout->thisMaskScan = 0; 278 tmpReadout->this WeightScan = 0;282 tmpReadout->thisVarianceScan = 0; 279 283 280 284 tmpReadout->lastImageScan = 0; 281 285 tmpReadout->lastMaskScan = 0; 282 tmpReadout->last WeightScan = 0;286 tmpReadout->lastVarianceScan = 0; 283 287 284 288 tmpReadout->forceScan = false; -
branches/pap_branch_20090128/psModules/src/camera/pmFPA.h
r19013 r21211 6 6 * @author Eugene Magnier, IfA 7 7 * 8 * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 8-08-12 02:51:20$8 * @version $Revision: 1.25.26.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii 11 11 */ … … 113 113 /// 114 114 /// A readout corresponds to an individual read of a cell (e.g., a single image as part of a video sequence, 115 /// or one of multiple coadds). It contains the actual pixels used in analysis (along with mask and weight115 /// or one of multiple coadds). It contains the actual pixels used in analysis (along with mask and variance 116 116 /// maps). When reading from a FITS file, the images are subimages (from CELL.TRIMSEC) of the pixels read 117 117 /// from the appropriate HDU (at the FPA, chip or cell level). The readout also contains a list of bias … … 123 123 psImage *image; ///< Imaging area of readout (corresponds to CELL.TRIMSEC region) 124 124 psImage *mask; ///< Mask of input image (corresponds to CELL.TRIMSEC region) 125 psImage *weight; ///< Weight of input image (corresponds to CELL.TRIMSEC region) 125 psImage *variance; ///< Variance of input image (corresponds to CELL.TRIMSEC region) 126 psKernel *covariance; ///< Covariance pseudo-matrix (covariance factors for single pixel) 126 127 psList *bias; ///< List of bias (prescan/overscan) images 127 128 psMetadata *analysis; ///< Readout-level analysis metadata … … 130 131 bool file_exists; ///< Does the file for this readout exist (read case only)? 131 132 bool data_exists; ///< Does the data for this readout exist (read case only)? 132 int thisImageScan; ///< start scan for next/current read 133 int lastImageScan; ///< start scan of the last read 134 int thisMaskScan; ///< start scan for next/current read 135 int lastMaskScan; ///< start scan of the last read 136 int this WeightScan; ///< start scan for next/current read137 int last WeightScan; ///< start scan of the last read133 int thisImageScan; ///< start scan for next/current read of image 134 int lastImageScan; ///< start scan of the last read of image 135 int thisMaskScan; ///< start scan for next/current read of mask 136 int lastMaskScan; ///< start scan of the last read of mask 137 int thisVarianceScan; ///< start scan for next/current read of variance 138 int lastVarianceScan; ///< start scan of the last read of variance 138 139 bool forceScan; ///< Force pmFPARead to obey the above commanded this and last scans. 139 140 } pmReadout; … … 217 218 } \ 218 219 } \ 219 psImage *weight = (READOUT)->weight; /* Weight map pixels */ \ 220 if (weight) { \ 221 PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL); \ 222 if ((numCols != 0 || numRows != 0) && (weight->numCols != numCols || weight->numRows != numRows)) { \ 223 psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Weight in readout %s has wrong size (%dx%d vs %dx%d)", \ 224 #READOUT, weight->numCols, weight->numRows, numCols, numRows); \ 220 psImage *variance = (READOUT)->variance; /* Variance map pixels */ \ 221 if (variance) { \ 222 PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL); \ 223 if ((numCols != 0 || numRows != 0) && \ 224 (variance->numCols != numCols || variance->numRows != numRows)) { \ 225 psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ 226 "Variance in readout %s has wrong size (%dx%d vs %dx%d)", \ 227 #READOUT, variance->numCols, variance->numRows, numCols, numRows); \ 225 228 return RETVAL; \ 226 229 } \ … … 236 239 PS_ASSERT_IMAGE_NON_NULL((READOUT)->mask, RETVAL); 237 240 238 /// Assert that a readout contains a weight map 239 #define PM_ASSERT_READOUT_WEIGHT(READOUT, RETVAL) \ 240 PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL); 241 /// Assert that a readout contains a variance map 242 #define PM_ASSERT_READOUT_VARIANCE(READOUT, RETVAL) \ 243 PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL); 244 245 /// Assert that a readout contains a covariance matrix 246 #define PM_ASSERT_READOUT_COVARIANCE(READOUT, RETVAL) \ 247 PS_ASSERT_KERNEL_NON_NULL((READOUT)->covariance, RETVAL); 241 248 242 249 /// @} -
branches/pap_branch_20090128/psModules/src/camera/pmFPAConstruct.c
r17911 r21211 1499 1499 psImage *image = readout->image; // The image 1500 1500 psImage *mask = readout->mask; // The mask 1501 psImage * weight = readout->weight; // The weight1501 psImage *variance = readout->variance; // The variance 1502 1502 psList *bias = readout->bias; // The list of bias images 1503 1503 if (image) { 1504 1504 INDENT(fd, 7); 1505 fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n", image->col0, image->col0 + 1506 image->numCols, image->row0, image->row0 + image->numRows, image->numCols, 1507 image->numRows); 1505 fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n", 1506 image->col0, image->col0 + image->numCols, 1507 image->row0, image->row0 + image->numRows, 1508 image->numCols, image->numRows); 1508 1509 } 1509 1510 if (bias) { … … 1512 1513 while ((biasImage = psListGetAndIncrement(biasIter))) { 1513 1514 INDENT(fd, 7); 1514 fprintf(fd, "Bias: [%d:%d,%d:%d] (%dx%d)\n", biasImage->col0, 1515 biasImage->col0 + biasImage->numCols, biasImage->row0, 1516 biasImage->row0 + biasImage->numRows, biasImage->numCols, biasImage->numRows); 1515 fprintf(fd, "Bias: [%d:%d,%d:%d] (%dx%d)\n", 1516 biasImage->col0, biasImage->col0 + biasImage->numCols, 1517 biasImage->row0, biasImage->row0 + biasImage->numRows, 1518 biasImage->numCols, biasImage->numRows); 1517 1519 } 1518 1520 psFree(biasIter); … … 1520 1522 if (mask) { 1521 1523 INDENT(fd, 7); 1522 fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n", mask->col0, mask->col0 + 1523 mask->numCols, mask->row0, mask->row0 + mask->numRows, mask->numCols, 1524 mask->numRows); 1524 fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n", 1525 mask->col0, mask->col0 + mask->numCols, 1526 mask->row0, mask->row0 + mask->numRows, 1527 mask->numCols, mask->numRows); 1525 1528 } 1526 if ( weight) {1529 if (variance) { 1527 1530 INDENT(fd, 7); 1528 fprintf(fd, "Weight: [%d:%d,%d:%d] (%dx%d)\n", weight->col0, weight->col0 + 1529 weight->numCols, weight->row0, weight->row0 + weight->numRows, weight->numCols, 1530 weight->numRows); 1531 fprintf(fd, "Variance: [%d:%d,%d:%d] (%dx%d)\n", 1532 variance->col0, variance->col0 + variance->numCols, 1533 variance->row0, variance->row0 + variance->numRows, 1534 variance->numCols, variance->numRows); 1531 1535 } 1532 1536 } // Iterating over cell -
branches/pap_branch_20090128/psModules/src/camera/pmFPACopy.c
r20634 r21211 224 224 targetReadout->data_exists = sourceReadout->data_exists; 225 225 226 // Copy all three image components (image, mask, weight) 227 readoutCopyComponent(&targetReadout->image, sourceReadout->image, binning, xFlip, yFlip, pixels); 228 readoutCopyComponent(&targetReadout->mask, sourceReadout->mask, binning, xFlip, yFlip, pixels); 229 readoutCopyComponent(&targetReadout->weight, sourceReadout->weight, binning, xFlip, yFlip, pixels); 226 // Copy all three image components (image, mask, variance) 227 readoutCopyComponent(&targetReadout->image, sourceReadout->image, binning, xFlip, yFlip, pixels); 228 readoutCopyComponent(&targetReadout->mask, sourceReadout->mask, binning, xFlip, yFlip, pixels); 229 readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip, 230 pixels); 230 231 231 232 // Copy bias … … 328 329 int xSize = psMetadataLookupS32(NULL, source->concepts, "CELL.XSIZE"); // CELL.XSIZE of source 329 330 330 psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source");331 psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source"); 331 332 332 333 if (sourceBin == 0) { … … 478 479 pmCell *targetCell = pmCellAlloc (targetChip, cellName); 479 480 int xParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.XPARITY"); // Target x parity 480 psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target");481 psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target"); 481 482 psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "", xParityTarget); 482 483 int yParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.YPARITY"); // Target y parity 483 psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target");484 psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target"); 484 485 psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "", yParityTarget); 485 486 if (!cellCopy(targetCell, sourceCell, true, 1, 1)) { -
branches/pap_branch_20090128/psModules/src/camera/pmFPAExtent.c
r20635 r21211 16 16 } 17 17 if (!image) { 18 image = readout-> weight;18 image = readout->variance; 19 19 } 20 20 … … 23 23 24 24 if (!image) { 25 pmHDU *hdu = pmHDUFromReadout (readout);26 if (hdu && hdu->header) {27 bool status;28 xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1");29 if (!status) {30 xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1");31 if (!status) return NULL;32 } 33 ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2");34 if (!status) {35 ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2");36 if (!status) return NULL;37 } 38 } else {25 pmHDU *hdu = pmHDUFromReadout (readout); 26 if (hdu && hdu->header) { 27 bool status; 28 xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1"); 29 if (!status) { 30 xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1"); 31 if (!status) return NULL; 32 } 33 ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2"); 34 if (!status) { 35 ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2"); 36 if (!status) return NULL; 37 } 38 } else { 39 39 // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent) 40 xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE");41 ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE");42 }40 xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE"); 41 ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE"); 42 } 43 43 return psRegionAlloc(0, xSize, 0, ySize); 44 44 } … … 70 70 // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent) 71 71 if (readouts->n == 0) { 72 int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");73 int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");74 cellExtent->x0 = 0;75 cellExtent->x1 = xSize;76 cellExtent->y0 = 0;77 cellExtent->y1 = ySize;72 int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE"); 73 int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE"); 74 cellExtent->x0 = 0; 75 cellExtent->x1 = xSize; 76 cellExtent->y0 = 0; 77 cellExtent->y1 = ySize; 78 78 } 79 79 … … 97 97 // CELL.X0,Y0 are the coordinate of the amp on the chip, subtract size if parity flipped 98 98 if (xParityCell > 0) { 99 cellExtent->x0 += x0;100 cellExtent->x1 += x0;99 cellExtent->x0 += x0; 100 cellExtent->x1 += x0; 101 101 } else { 102 float x0Cell = x0 - cellExtent->x1;103 float x1Cell = x0 - cellExtent->x0;104 cellExtent->x0 = x0Cell;105 cellExtent->x1 = x1Cell;102 float x0Cell = x0 - cellExtent->x1; 103 float x1Cell = x0 - cellExtent->x0; 104 cellExtent->x0 = x0Cell; 105 cellExtent->x1 = x1Cell; 106 106 } 107 107 if (yParityCell > 0) { 108 cellExtent->y0 += y0;109 cellExtent->y1 += y0;108 cellExtent->y0 += y0; 109 cellExtent->y1 += y0; 110 110 } else { 111 float y0Cell = y0 - cellExtent->y1;112 float y1Cell = y0 - cellExtent->y0;113 cellExtent->y0 = y0Cell;114 cellExtent->y1 = y1Cell;111 float y0Cell = y0 - cellExtent->y1; 112 float y1Cell = y0 - cellExtent->y0; 113 cellExtent->y0 = y0Cell; 114 cellExtent->y1 = y1Cell; 115 115 } 116 116 -
branches/pap_branch_20090128/psModules/src/camera/pmFPAMaskWeight.c
r21183 r21211 49 49 } 50 50 51 // Create the parent weightimages that reside in the HDU52 static void createParent Weights(pmHDU *hdu // The HDU for which to create51 // Create the parent variance images that reside in the HDU 52 static void createParentVariances(pmHDU *hdu // The HDU for which to create 53 53 ) 54 54 { … … 58 58 // Generate the parent mask images 59 59 psArray *images = hdu->images; // Array of images 60 psArray * weights = hdu->weights; // Array of weightimages61 if (! weights) {62 weights = psArrayAlloc(images->n);63 hdu-> weights = weights;60 psArray *variances = hdu->variances; // Array of variance images 61 if (!variances) { 62 variances = psArrayAlloc(images->n); 63 hdu->variances = variances; 64 64 } 65 65 66 66 for (long i = 0; i < images->n; i++) { 67 67 psImage *image = images->data[i]; // The image for this readout 68 if (!image || weights->data[i]) {68 if (!image || variances->data[i]) { 69 69 continue; 70 70 } 71 weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);71 variances->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); 72 72 } 73 73 … … 199 199 } 200 200 201 bool pmReadoutSet Weight(pmReadout *readout, bool poisson)201 bool pmReadoutSetVariance(pmReadout *readout, bool poisson) 202 202 { 203 203 PS_ASSERT_PTR_NON_NULL(readout, false); … … 209 209 float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain 210 210 if (!mdok || isnan(gain)) { 211 psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set weight.\n");211 psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n"); 212 212 return false; 213 213 } 214 214 float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise 215 215 if (!mdok || isnan(readnoise)) { 216 psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set weight.\n");216 psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n"); 217 217 return false; 218 218 } … … 223 223 224 224 if (poisson) { 225 // Set weightimage to the variance in ADU = f/g + rn^2225 // Set variance image to the variance in ADU = f/g + rn^2 226 226 psImage *image = readout->image; // The image pixels 227 readout-> weight = (psImage*)psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32));228 229 // a negative weight is non-sensical. if the image value drops below 1, the weightmust be 1.230 readout-> weight = (psImage*)psUnaryOp(readout->weight, readout->weight, "abs");231 readout-> weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "max",227 readout->variance = (psImage*)psBinaryOp(readout->variance, image, "/", psScalarAlloc(gain, PS_TYPE_F32)); 228 229 // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1. 230 readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs"); 231 readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max", 232 232 psScalarAlloc(1, PS_TYPE_F32)); 233 233 } else { 234 234 // Just use the read noise 235 if (!readout-> weight) {236 readout-> weight= psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);237 } 238 psImageInit(readout-> weight, 0.0);239 } 240 241 readout-> weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "+",235 if (!readout->variance) { 236 readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); 237 } 238 psImageInit(readout->variance, 0.0); 239 } 240 241 readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", 242 242 psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); 243 243 … … 245 245 } 246 246 247 // this function creates the weight pixels, or uses the existing weightpixels. it will set248 // the noise pixel values only if the weightimage is not supplied249 bool pmReadoutGenerate Weight(pmReadout *readout, bool poisson)247 // this function creates the variance pixels, or uses the existing variance pixels. it will set 248 // the noise pixel values only if the variance image is not supplied 249 bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson) 250 250 { 251 251 PS_ASSERT_PTR_NON_NULL(readout, false); … … 254 254 bool mdok = true; // Status of MD lookup 255 255 256 // Create the weightimage if required257 if (readout-> weight)256 // Create the variance image if required 257 if (readout->variance) 258 258 return true; 259 259 260 260 psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section 261 261 if (!mdok || psRegionIsNaN(*trimsec)) { 262 psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set weight.\n");262 psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n"); 263 263 return false; 264 264 } … … 271 271 } 272 272 273 createParent Weights(hdu);273 createParentVariances(hdu); 274 274 275 275 // Need to identify which readout we're working with.... … … 280 280 } 281 281 282 psImage * weight = psImageSubset(hdu->weights->data[index], *trimsec); // The weightpixels283 if (! weight) {282 psImage *variance = psImageSubset(hdu->variances->data[index], *trimsec); // The variance pixels 283 if (!variance) { 284 284 psString trimsecString = psRegionToString(*trimsec); 285 psError(PS_ERR_UNKNOWN, false, "Unable to set weightfrom HDU with trimsec: %s.\n",285 psError(PS_ERR_UNKNOWN, false, "Unable to set variance from HDU with trimsec: %s.\n", 286 286 trimsecString); 287 287 psFree(trimsecString); 288 288 return false; 289 289 } 290 psImageInit( weight, 0);291 readout-> weight = weight;292 293 return pmReadoutSet Weight(readout, poisson);294 } 295 296 bool pmReadoutGenerateMask Weight(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)290 psImageInit(variance, 0); 291 readout->variance = variance; 292 293 return pmReadoutSetVariance(readout, poisson); 294 } 295 296 bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson) 297 297 { 298 298 PS_ASSERT_PTR_NON_NULL(readout, false); … … 301 301 302 302 success &= pmReadoutGenerateMask(readout, satMask, badMask); 303 success &= pmReadoutGenerate Weight(readout, poisson);303 success &= pmReadoutGenerateVariance(readout, poisson); 304 304 305 305 return success; 306 306 } 307 307 308 bool pmCellGenerateMask Weight(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)308 bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson) 309 309 { 310 310 PS_ASSERT_PTR_NON_NULL(cell, false); … … 314 314 for (int i = 0; i < readouts->n; i++) { 315 315 pmReadout *readout = readouts->data[i]; // The readout 316 success &= pmReadoutGenerateMask Weight(readout, poisson, satMask, badMask);316 success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask); 317 317 } 318 318 … … 321 321 322 322 323 bool pmReadout WeightRenormPixels(const pmReadout *readout, psImageMaskType maskVal,323 bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal, 324 324 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) 325 325 { 326 326 PM_ASSERT_READOUT_NON_NULL(readout, false); 327 327 PM_ASSERT_READOUT_IMAGE(readout, false); 328 PM_ASSERT_READOUT_ WEIGHT(readout, false);329 330 psImage *image = readout->image, *mask = readout->mask, * weight = readout->weight; // Readout parts328 PM_ASSERT_READOUT_VARIANCE(readout, false); 329 330 psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts 331 331 332 332 if (!psMemIncrRefCounter(rng)) { … … 334 334 } 335 335 336 psStats * weightStats = psStatsAlloc(meanStat);// Statistics for mean337 if (!psImageBackground( weightStats, NULL, weight, mask, maskVal, rng)) {338 psError(PS_ERR_UNKNOWN, false, "Unable to measure mean weightfor image");339 psFree( weightStats);336 psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean 337 if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) { 338 psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image"); 339 psFree(varianceStats); 340 340 psFree(rng); 341 341 return false; 342 342 } 343 float meanVariance = weightStats->robustMedian; // Mean variance344 psFree( weightStats);343 float meanVariance = varianceStats->robustMedian; // Mean variance 344 psFree(varianceStats); 345 345 346 346 psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean … … 356 356 357 357 float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be 358 psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weightmap by %f", correction);359 psBinaryOp( weight, weight, "*", psScalarAlloc(correction, PS_TYPE_F32));358 psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction); 359 psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32)); 360 360 361 361 return true; … … 363 363 364 364 365 bool pmReadout WeightRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,365 bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width, 366 366 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) 367 367 { 368 368 PM_ASSERT_READOUT_NON_NULL(readout, false); 369 369 PM_ASSERT_READOUT_IMAGE(readout, false); 370 PM_ASSERT_READOUT_ WEIGHT(readout, false);370 PM_ASSERT_READOUT_VARIANCE(readout, false); 371 371 372 372 if (!psMemIncrRefCounter(rng)) { … … 374 374 } 375 375 376 psImage *image = readout->image, *mask = readout->mask, * weight = readout->weight; // Readout images376 psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images 377 377 378 378 // Measure background … … 424 424 float sumNoise = 0.0; // Sum for noise measurement 425 425 float sumSource = 0.0; // Sum for source measurement 426 float sum Weight = 0.0; // Sum for weightmeasurement426 float sumVariance = 0.0; // Sum for variance measurement 427 427 float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels 428 428 for (int v = 0, y = yPix - size; v < fullSize; v++, y++) { 429 429 float xSumNoise = 0.0; // Sum for noise measurement in x 430 430 float xSumSource = 0.0; // Sum for source measurement in x 431 float xSum Weight = 0.0; // Sum for weightmeasurement in x431 float xSumVariance = 0.0; // Sum for variance measurement in x 432 432 float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x 433 433 float yGauss = gauss->data.F32[v]; // Value of Gaussian in y … … 441 441 xSumNoise += value * xGauss; 442 442 xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss; 443 xSum Weight += weight->data.F32[y][x] * xGauss2;443 xSumVariance += variance->data.F32[y][x] * xGauss2; 444 444 xSumGauss += xGauss; 445 445 xSumGauss2 += xGauss2; … … 448 448 sumNoise += xSumNoise * yGauss; 449 449 sumSource += xSumSource * yGauss; 450 sum Weight += xSumWeight* yGauss2;450 sumVariance += xSumVariance * yGauss2; 451 451 sumGauss += xSumGauss * yGauss; 452 452 sumGauss2 += xSumGauss2 * yGauss2; … … 454 454 455 455 photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) && 456 isfinite(sum Weight) && sumGauss > 0 && sumGauss2 > 0) ?456 isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ? 457 457 0 : 0xFF); 458 458 459 459 float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise 460 460 float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source 461 float smooth Weight = sumWeight / sumGauss2; // Value of smoothed weightpixel461 float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel 462 462 463 463 noise->data.F32[i] = smoothImageNoise; 464 464 source->data.F32[i] = smoothImageSource; 465 guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smooth Weight: 0.0;465 guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0; 466 466 psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n", 467 i, xPix, yPix, smoothImageNoise, smoothImageSource, smooth Weight);467 i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance); 468 468 } 469 469 psFree(gauss); … … 516 516 psFree(photMask); 517 517 518 psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weightmap by %f", meanRatio);519 psBinaryOp( weight, weight, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));518 psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio); 519 psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32)); 520 520 521 521 return true; … … 523 523 524 524 525 bool pmReadout WeightRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,525 bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat, 526 526 psStatsOptions stdevStat, int width, psRandom *rng) 527 527 { 528 528 PM_ASSERT_READOUT_NON_NULL(readout, false); 529 529 PM_ASSERT_READOUT_IMAGE(readout, false); 530 PM_ASSERT_READOUT_ WEIGHT(readout, false);530 PM_ASSERT_READOUT_VARIANCE(readout, false); 531 531 PS_ASSERT_INT_POSITIVE(width, false); 532 532 … … 535 535 } 536 536 537 psImage *image = readout->image, *mask = readout->mask, * weight = readout->weight; // Readout images537 psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images 538 538 int numCols = image->numCols, numRows = image->numRows; // Size of images 539 539 int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions … … 554 554 psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest 555 555 psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels 556 psImage *sub Weight = psImageSubset(weight, region); // Sub image of the weightpixels556 psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels 557 557 psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels 558 558 559 559 if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) || 560 !psImageBackground(meanStats, &buffer, sub Weight, subMask, maskVal, rng)) {560 !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) { 561 561 // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it. 562 562 psString regionStr = psRegionToString(region); // String with region … … 564 564 psFree(regionStr); 565 565 psErrorClear(); 566 psImageInit(sub Weight, NAN);566 psImageInit(subVariance, NAN); 567 567 if (subMask) { 568 568 psImageInit(subMask, maskVal); … … 574 574 "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n", 575 575 xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar); 576 psBinaryOp(sub Weight, subWeight, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));576 psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32)); 577 577 } 578 578 579 579 psFree(subImage); 580 psFree(sub Weight);580 psFree(subVariance); 581 581 psFree(subMask); 582 582 } … … 597 597 598 598 psImage *image = readout->image; // Readout's image 599 psImage * weight = readout->weight; // Readout's weight599 psImage *variance = readout->variance; // Readout's variance 600 600 int numCols = image->numCols, numRows = image->numRows; // Size of image 601 601 … … 607 607 for (int y = 0; y < numRows; y++) { 608 608 for (int x = 0; x < numCols; x++) { 609 if (!isfinite(image->data.F32[y][x]) || ( weight && !isfinite(weight->data.F32[y][x]))) {609 if (!isfinite(image->data.F32[y][x]) || (variance && !isfinite(variance->data.F32[y][x]))) { 610 610 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal; 611 611 } … … 627 627 psImageMaskType **maskData = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Dereference mask 628 628 psF32 **imageData = readout->image->data.F32;// Dereference image 629 psF32 ** weightData = readout->weight ? readout->weight->data.F32 : NULL; // Dereference weightmap629 psF32 **varianceData = readout->variance ? readout->variance->data.F32 : NULL; // Dereference variance map 630 630 631 631 for (int y = 0; y < numRows; y++) { … … 633 633 if (maskData[y][x] & maskVal) { 634 634 imageData[y][x] = NAN; 635 if ( weightData) {636 weightData[y][x] = NAN;635 if (varianceData) { 636 varianceData[y][x] = NAN; 637 637 } 638 638 } … … 656 656 psImage *image = readout->image; // Image 657 657 psImage *mask = readout->mask; // Mask 658 psImage * weight = readout->weight; // Weightmap659 660 psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal,658 psImage *variance = readout->variance; // Variance map 659 660 psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal, 661 661 NAN, NAN, maskBad, maskPoor, poorFrac, 0); 662 662 interp->shifting = false; // Turn off "exact shifts" so we get proper interpolation … … 666 666 psPixels *pixels = psPixelsAllocEmpty(PIXELS_BUFFER); // Pixels that have been interpolated 667 667 psVector *imagePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for image 668 psVector * weightPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for weight668 psVector *variancePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for variance 669 669 psVector *maskPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_IMAGE_MASK); // Corresponding values for mask 670 670 // NOTE: maskPix carries the actual image mask values -- do NOT use … … 675 675 for (int x = 0; x < numCols; x++) { 676 676 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { 677 double imageValue, weightValue; // Image and weightvalue from interpolation677 double imageValue, varianceValue; // Image and variance value from interpolation 678 678 psImageMaskType maskValue = 0; // Mask value from interpolation 679 psImageInterpolateStatus status = psImageInterpolate(&imageValue, & weightValue, &maskValue, x, y, interp);679 psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp); 680 680 if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) { 681 681 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y); … … 683 683 psFree(pixels); 684 684 psFree(imagePix); 685 psFree( weightPix);685 psFree(variancePix); 686 686 psFree(maskPix); 687 687 return false; … … 694 694 pixels = psPixelsAdd(pixels, PIXELS_BUFFER, x, y); 695 695 imagePix = psVectorExtend(imagePix, PIXELS_BUFFER, 1); 696 weightPix = psVectorExtend(weightPix, PIXELS_BUFFER, 1);696 variancePix = psVectorExtend(variancePix, PIXELS_BUFFER, 1); 697 697 maskPix = psVectorExtend(maskPix, PIXELS_BUFFER, 1); 698 698 imagePix->data.F32[numBad] = imageValue; 699 weightPix->data.F32[numBad] = weightValue;699 variancePix->data.F32[numBad] = varianceValue; 700 700 maskPix->data.PS_TYPE_IMAGE_MASK_DATA[numBad] = maskValue; 701 701 numBad++; … … 709 709 int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of pixel 710 710 image->data.F32[y][x] = imagePix->data.F32[i]; 711 weight->data.F32[y][x] = weightPix->data.F32[i];711 variance->data.F32[y][x] = variancePix->data.F32[i]; 712 712 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskPix->data.PS_TYPE_IMAGE_MASK_DATA[i]; 713 713 } … … 715 715 psFree(pixels); 716 716 psFree(imagePix); 717 psFree( weightPix);717 psFree(variancePix); 718 718 psFree(maskPix); 719 719 -
branches/pap_branch_20090128/psModules/src/camera/pmFPAMaskWeight.h
r21183 r21211 5 5 * @author Eugene Magnier, IfA 6 6 * 7 * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $8 * @date $Date: 2009-01-2 7 06:39:38$7 * @version $Revision: 1.17.2.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2009-01-29 00:33:51 $ 9 9 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii 10 10 */ … … 47 47 48 48 49 /// Set a temporary readout weightmap using CELL.GAIN and CELL.READNOISE49 /// Set a temporary readout variance map using CELL.GAIN and CELL.READNOISE 50 50 /// 51 /// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain52 /// (CELL. GAIN) and read noise (CELL.READNOISE). The weight map that is produced within the readout is53 /// t emporary --- it is not added to the HDU. This is intended for when the user is iterating using54 /// pmReadoutReadNext, in which case the HDUcan't be generated.55 bool pmReadoutSet Weight(pmReadout *readout, ///< Readout for which to set weight56 bool poisson ///< Use poisson weights(in addition to read noise)?57 );51 /// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise 52 /// (CELL.READNOISE). The weight map that is produced within the readout is temporary --- it is not added to 53 /// the HDU. This is intended for when the user is iterating using pmReadoutReadNext, in which case the HDU 54 /// can't be generated. 55 bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance 56 bool poisson ///< Include poisson variance (in addition to read noise)? 57 ); 58 58 59 59 /// Generate a readout mask (suitable for output) using CELL.SATURATION and CELL.BAD … … 66 66 ); 67 67 68 /// Generate a weightmap (suitable for output) using CELL.GAIN and CELL.READNOISE68 /// Generate a variance map (suitable for output) using CELL.GAIN and CELL.READNOISE 69 69 /// 70 /// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain71 /// (CELL. GAIN) and read noise (CELL.READNOISE). The weight map that is produced within the readout is72 /// suitable for output (completewith HDU entry). This is intended for most operations.73 bool pmReadoutGenerate Weight(pmReadout *readout, ///< Readout for which to generate weight74 bool poisson ///< Use poisson weights(in addition to read noise)?75 );70 /// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise 71 /// (CELL.READNOISE). The variance map that is produced within the readout is suitable for output (complete 72 /// with HDU entry). This is intended for most operations. 73 bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance 74 bool poisson ///< Include poisson variance (in addition to read noise)? 75 ); 76 76 77 /// Generate mask and weightmap for a readout77 /// Generate mask and variance map for a readout 78 78 /// 79 /// Calls pmReadoutGenerateMask and pmReadoutGenerate Weightfor the readout80 bool pmReadoutGenerateMask Weight(pmReadout *readout, ///< Readout for which to generate mask and weights81 psImageMaskType sat, ///< Mask value to give saturated pixels82 psImageMaskType bad, ///< Mask value to give bad (low) pixels83 bool poisson ///< Use poisson weights(in addition to read noise)?84 );79 /// Calls pmReadoutGenerateMask and pmReadoutGenerateVariance for the readout 80 bool pmReadoutGenerateMaskVariance(pmReadout *readout, ///< Readout for which to generate mask and variance 81 psImageMaskType sat, ///< Mask value to give saturated pixels 82 psImageMaskType bad, ///< Mask value to give bad (low) pixels 83 bool poisson ///< Include poisson variance (in addition to read noise)? 84 ); 85 85 86 /// Generate mask and weightmaps for all readouts within a cell86 /// Generate mask and variance maps for all readouts within a cell 87 87 /// 88 /// Calls pmReadoutGenerateMask Weightfor each readout within the cell.89 bool pmCellGenerateMask Weight(pmCell *cell, ///< Cell for which to generate mask and weights90 psImageMaskType sat, ///< Mask value to give saturated pixels91 psImageMaskType bad, ///< Mask value to give bad (low) pixels92 bool poisson ///< Use poisson weights(in addition to read noise)?93 );88 /// Calls pmReadoutGenerateMaskVariance for each readout within the cell. 89 bool pmCellGenerateMaskVariance(pmCell *cell, ///< Cell for which to generate mask and variance 90 psImageMaskType sat, ///< Mask value to give saturated pixels 91 psImageMaskType bad, ///< Mask value to give bad (low) pixels 92 bool poisson ///< Include poisson variance (in addition to read noise)? 93 ); 94 94 95 /// Renormalise the weightmap to match the actual pixel variance95 /// Renormalise the variance map to match the actual pixel variance 96 96 /// 97 /// The weight (variance)map is adjusted so that the mean matches the actual pixel variance in the image98 bool pmReadout WeightRenormPixels(97 /// The variance map is adjusted so that the mean matches the actual pixel variance in the image 98 bool pmReadoutVarianceRenormPixels( 99 99 const pmReadout *readout, ///< Readout to normalise 100 100 psImageMaskType maskVal, ///< Value to mask … … 104 104 ); 105 105 106 /// Renormalise the weightmap to match the actual photometry variance106 /// Renormalise the variance map to match the actual photometry variance 107 107 /// 108 /// The weight (variance)map is adjusted so that the actual significance of fake sources matches the108 /// The variance map is adjusted so that the actual significance of fake sources matches the 109 109 /// guestimated significance 110 bool pmReadout WeightRenormPhot(110 bool pmReadoutVarianceRenormPhot( 111 111 const pmReadout *readout, ///< Readout to normalise 112 psImageMaskType maskVal, ///< Value to mask112 psImageMaskType maskVal, ///< Value to mask 113 113 int num, ///< Number of instances to measure over the image 114 114 float width, ///< Photometry width … … 118 118 ); 119 119 120 /// Renormalise the weightmap to match the actual variance120 /// Renormalise the variance map to match the actual variance 121 121 /// 122 122 /// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for 123 123 /// that patch corresponds. 124 bool pmReadout WeightRenorm(const pmReadout *readout, // Readout to normalise125 psImageMaskType maskVal, // Value to mask126 psStatsOptions meanStat, // Statistic to measure the mean (of the variance map)127 psStatsOptions stdevStat, // Statistic to measure the stdev (of the image)128 int width, // Width of patch (pixels)129 psRandom *rng // Random number generator (for sub-sampling images)124 bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise 125 psImageMaskType maskVal, // Value to mask 126 psStatsOptions meanStat, // Statistic to measure the mean (of the variance map) 127 psStatsOptions stdevStat, // Statistic to measure the stdev (of the image) 128 int width, // Width of patch (pixels) 129 psRandom *rng // Random number generator (for sub-sampling images) 130 130 ); 131 131 … … 133 133 /// 134 134 /// Since unmasked non-finite pixels can occur (e.g., by out-of-range in quantisation), it is sometimes 135 /// necessary to mask them explicitly. Non-finite pixels in the image or weighthave their mask OR-ed with135 /// necessary to mask them explicitly. Non-finite pixels in the image or variance have their mask OR-ed with 136 136 /// the provided value. 137 137 bool pmReadoutMaskNonfinite(pmReadout *readout, ///< Readout to mask … … 139 139 ); 140 140 141 /// Apply a mask to the image and weightmap141 /// Apply a mask to the image and variance map 142 142 /// 143 143 /// Unfortunately, image subtraction may result in a bi-modal image in masked areas, which can upset image 144 144 /// statistics (very important for quantising images so that a product can be written out!). This function 145 /// sets masked areas to NAN in the image and weight.145 /// sets masked areas to NAN in the image and variance. 146 146 bool pmReadoutMaskApply(pmReadout *readout, ///< Readout to mask 147 147 psImageMaskType maskVal ///< Mask value for which to apply mask -
branches/pap_branch_20090128/psModules/src/camera/pmFPAMosaic.c
r21183 r21211 554 554 // We have to do all of the hard work ourselves 555 555 switch (type) { 556 FILL_IN(U8);557 FILL_IN(U16);558 FILL_IN(U32);559 FILL_IN(U64);560 FILL_IN(S8);561 FILL_IN(S16);562 FILL_IN(S32);563 FILL_IN(S64);564 FILL_IN(F32);565 FILL_IN(F64);556 FILL_IN(U8); 557 FILL_IN(U16); 558 FILL_IN(U32); 559 FILL_IN(U64); 560 FILL_IN(S8); 561 FILL_IN(S16); 562 FILL_IN(S32); 563 FILL_IN(S64); 564 FILL_IN(F32); 565 FILL_IN(F64); 566 566 default: 567 567 psAbort("Should never get here.\n"); … … 575 575 static bool addCell(psArray *images, // Array of images 576 576 psArray *masks, // Array of masks 577 psArray * weights, // Array of weights577 psArray *variances, // Array of variances 578 578 psVector *x0, // Array of X0 579 579 psVector *y0, // Array of Y0 … … 604 604 images = psArrayRealloc(images, index + CELL_LIST_BUFFER); 605 605 masks = psArrayRealloc(masks, index + CELL_LIST_BUFFER); 606 weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER);606 variances = psArrayRealloc(variances, index + CELL_LIST_BUFFER); 607 607 x0 = psVectorRealloc(x0, index+ CELL_LIST_BUFFER); 608 608 y0 = psVectorRealloc(y0, index+ CELL_LIST_BUFFER); … … 615 615 images->n = index + 1; 616 616 masks->n = index + 1; 617 weights->n = index + 1;617 variances->n = index + 1; 618 618 x0->n = index + 1; 619 619 y0->n = index + 1; … … 726 726 // The images to put into the mosaic 727 727 images->data[index] = psMemIncrRefCounter(readout->image); 728 weights->data[index] = psMemIncrRefCounter(readout->weight);728 variances->data[index] = psMemIncrRefCounter(readout->variance); 729 729 masks->data[index] = psMemIncrRefCounter(readout->mask); 730 730 … … 740 740 static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned 741 741 psImage **mosaicMask, // The mosaic mask, to be returned 742 psImage **mosaic Weights, // The mosaic weights, to be returned742 psImage **mosaicVariances, // The mosaic variances, to be returned 743 743 int *xBinChip, int *yBinChip, // The binning in x and y, to be returned 744 744 const pmChip *chip, // Chip to mosaic … … 749 749 assert(mosaicImage); 750 750 assert(mosaicMask); 751 assert(mosaic Weights);751 assert(mosaicVariances); 752 752 assert(xBinChip); 753 753 assert(yBinChip); … … 756 756 757 757 psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked 758 psArray * weights = psArrayAlloc(0); // Array of weightimages to be mosaicked758 psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked 759 759 psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked 760 760 psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates … … 802 802 continue; 803 803 } 804 allGood &= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,804 allGood &= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip, 805 805 cell, xBinChip, yBinChip, false, x0Target, y0Target, 806 806 xParityTarget, yParityTarget); … … 826 826 if (allGood) { 827 827 *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); 828 *mosaic Weights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);828 *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); 829 829 *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank); 830 830 } … … 832 832 // Clean up 833 833 psFree(images); 834 psFree( weights);834 psFree(variances); 835 835 psFree(masks); 836 836 psFree(xFlip); … … 847 847 static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned 848 848 psImage **mosaicMask, // The mosaic mask, to be returned 849 psImage **mosaic Weights, // The mosaic weights, to be returned849 psImage **mosaicVariances, // The mosaic variances, to be returned 850 850 int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned 851 851 const pmFPA *fpa, // FPA to mosaic … … 857 857 assert(mosaicImage); 858 858 assert(mosaicMask); 859 assert(mosaic Weights);859 assert(mosaicVariances); 860 860 assert(xBinFPA); 861 861 assert(yBinFPA); … … 865 865 866 866 psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked 867 psArray * weights = psArrayAlloc(0); // Array of weightimages to be mosaicked867 psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked 868 868 psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked 869 869 psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates … … 941 941 continue; 942 942 } 943 allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,943 allGood |= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip, 944 944 cell, xBinFPA, yBinFPA, true, x0Target, y0Target, 945 945 xParityTarget, yParityTarget); … … 960 960 if (allGood) { 961 961 *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); 962 *mosaic Weights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);962 *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); 963 963 *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank); 964 964 } … … 966 966 // Clean up 967 967 psFree(images); 968 psFree( weights);968 psFree(variances); 969 969 psFree(masks); 970 970 psFree(xFlip); … … 1025 1025 psImage *mosaicImage = NULL; // The mosaic image 1026 1026 psImage *mosaicMask = NULL; // The mosaic mask 1027 psImage *mosaic Weights = NULL; // The mosaic weights1027 psImage *mosaicVariances = NULL; // The mosaic variances 1028 1028 1029 1029 // Find the HDU … … 1051 1051 } 1052 1052 } 1053 if (hdu-> weights) {1054 mosaic Weights = psImageSubset(hdu->weights->data[0], bounds);1055 if (!mosaic Weights) {1056 psError(PS_ERR_UNKNOWN, false, "Unable to select weightpixels.\n");1053 if (hdu->variances) { 1054 mosaicVariances = psImageSubset(hdu->variances->data[0], bounds); 1055 if (!mosaicVariances) { 1056 psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n"); 1057 1057 return false; 1058 1058 } … … 1061 1061 // Case 2 --- we need to mosaic by cut and paste 1062 1062 psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n"); 1063 if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaic Weights, &xBin, &yBin, source, targetCell, blank)) {1063 if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) { 1064 1064 psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n"); 1065 1065 return false; … … 1094 1094 newReadout->image = mosaicImage; 1095 1095 newReadout->mask = mosaicMask; 1096 newReadout-> weight = mosaicWeights;1096 newReadout->variance = mosaicVariances; 1097 1097 psFree(newReadout); // Drop reference 1098 1098 … … 1262 1262 psImage *mosaicImage = NULL; // The mosaic image 1263 1263 psImage *mosaicMask = NULL; // The mosaic mask 1264 psImage *mosaic Weights = NULL; // The mosaic weights1264 psImage *mosaicVariances = NULL; // The mosaic variances 1265 1265 1266 1266 // Find the HDU … … 1275 1275 mosaicMask = psImageSubset(hdu->masks->data[0], *fpaRegion); 1276 1276 } 1277 if (hdu-> weights) {1278 mosaic Weights = psImageSubset(hdu->weights->data[0], *fpaRegion);1277 if (hdu->variances) { 1278 mosaicVariances = psImageSubset(hdu->variances->data[0], *fpaRegion); 1279 1279 } 1280 1280 } else { 1281 1281 // Case 2 --- we need to mosaic by cut and paste 1282 1282 psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n"); 1283 if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaic Weights, &xBin, &yBin, source,1283 if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, 1284 1284 targetChip, targetCell, blank)) { 1285 1285 psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n"); … … 1338 1338 newReadout->image = mosaicImage; 1339 1339 newReadout->mask = mosaicMask; 1340 newReadout-> weight = mosaicWeights;1340 newReadout->variance = mosaicVariances; 1341 1341 psFree(newReadout); // Drop reference 1342 1342 -
branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c
r21183 r21211 28 28 FPA_READ_TYPE_IMAGE, // Read image 29 29 FPA_READ_TYPE_MASK, // Read mask 30 FPA_READ_TYPE_ WEIGHT, // Read weightmap30 FPA_READ_TYPE_VARIANCE, // Read variance map 31 31 FPA_READ_TYPE_HEADER // Read header 32 32 } fpaReadType; … … 54 54 case FPA_READ_TYPE_MASK: 55 55 return readout->thisMaskScan; 56 case FPA_READ_TYPE_ WEIGHT:57 return readout->this WeightScan;56 case FPA_READ_TYPE_VARIANCE: 57 return readout->thisVarianceScan; 58 58 default: 59 59 psAbort("Unknown read type: %x\n", type); … … 74 74 readout->thisMaskScan = thisScan; 75 75 return readout->lastMaskScan; 76 case FPA_READ_TYPE_ WEIGHT:77 readout->this WeightScan = thisScan;78 return readout->last WeightScan;76 case FPA_READ_TYPE_VARIANCE: 77 readout->thisVarianceScan = thisScan; 78 return readout->lastVarianceScan; 79 79 default: 80 80 psAbort("Unknown read type: %x\n", type); … … 93 93 case FPA_READ_TYPE_MASK: 94 94 return readout->lastMaskScan; 95 case FPA_READ_TYPE_ WEIGHT:96 return readout->last WeightScan;95 case FPA_READ_TYPE_VARIANCE: 96 return readout->lastVarianceScan; 97 97 default: 98 98 psAbort("Unknown read type: %x\n", type); … … 113 113 readout->lastMaskScan = lastScan; 114 114 return readout->lastMaskScan; 115 case FPA_READ_TYPE_ WEIGHT:116 readout->last WeightScan = lastScan;117 return readout->last WeightScan;115 case FPA_READ_TYPE_VARIANCE: 116 readout->lastVarianceScan = lastScan; 117 return readout->lastVarianceScan; 118 118 default: 119 119 psAbort("Unknown read type: %x\n", type); … … 132 132 case FPA_READ_TYPE_MASK: 133 133 return &readout->mask; 134 case FPA_READ_TYPE_ WEIGHT:135 return &readout-> weight;134 case FPA_READ_TYPE_VARIANCE: 135 return &readout->variance; 136 136 default: 137 137 psAbort("Unknown read type: %x\n", type); … … 350 350 *target = psImageSubset(image, region); 351 351 352 // Get the list of overscans: only for IMAGE types (no overscan for MASK and WEIGHT)352 // Get the list of overscans: only for IMAGE types (no overscan for MASK and VARIANCE) 353 353 if (type == FPA_READ_TYPE_IMAGE) { 354 354 if (readout->bias->n != 0) { … … 582 582 *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]); 583 583 584 // Read overscans only for "image" type --- weights and masks shouldn't record overscans584 // Read overscans only for "image" type --- variances and masks shouldn't record overscans 585 585 if (type == FPA_READ_TYPE_IMAGE) { 586 586 // Blow away existing data … … 608 608 } 609 609 610 // Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellRead Weight610 // Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadVariance 611 611 // Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts. 612 612 static bool cellRead(pmCell *cell, // Cell into which to read … … 640 640 dataPointer = hdu->masks; 641 641 break; 642 case FPA_READ_TYPE_ WEIGHT:643 hduReadFunc = pmHDURead Weight;644 dataPointer = hdu-> weights;642 case FPA_READ_TYPE_VARIANCE: 643 hduReadFunc = pmHDUReadVariance; 644 dataPointer = hdu->variances; 645 645 break; 646 646 default: … … 680 680 imageArray = hdu->masks; 681 681 break; 682 case FPA_READ_TYPE_ WEIGHT:683 imageArray = hdu-> weights;682 case FPA_READ_TYPE_VARIANCE: 683 imageArray = hdu->variances; 684 684 break; 685 685 default: … … 729 729 730 730 731 // Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipRead Weight731 // Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadVariance 732 732 // Iterates over component cells, reading each 733 733 static bool chipRead(pmChip *chip, // Chip into which to read … … 760 760 761 761 762 // Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPARead Weight762 // Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadVariance 763 763 // Iterates over component chips, reading each 764 764 static bool fpaRead(pmFPA *fpa, // FPA into which to read … … 1091 1091 1092 1092 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 1093 // Reading the weightmap1094 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 1095 1096 bool pmReadoutMore Weight(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)1093 // Reading the variance map 1094 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 1095 1096 bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) 1097 1097 { 1098 1098 PS_ASSERT_PTR_NON_NULL(readout, false); 1099 1099 PS_ASSERT_FITS_NON_NULL(fits, false); 1100 1100 1101 return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_ WEIGHT, config);1102 } 1103 1104 bool pmReadoutReadChunk Weight(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,1101 return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config); 1102 } 1103 1104 bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, 1105 1105 pmConfig *config) 1106 1106 { … … 1110 1110 PS_ASSERT_INT_NONNEGATIVE(numScans, false); 1111 1111 1112 return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_ WEIGHT, config);1113 } 1114 1115 bool pmReadoutRead Weight(pmReadout *readout, psFits *fits, int z, pmConfig *config)1112 return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config); 1113 } 1114 1115 bool pmReadoutReadVariance(pmReadout *readout, psFits *fits, int z, pmConfig *config) 1116 1116 { 1117 1117 PS_ASSERT_PTR_NON_NULL(readout, false); 1118 1118 PS_ASSERT_FITS_NON_NULL(fits, false); 1119 1119 1120 return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_ WEIGHT, config);1121 } 1122 1123 bool pmCellRead Weight(pmCell *cell, psFits *fits, pmConfig *config)1120 return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config); 1121 } 1122 1123 bool pmCellReadVariance(pmCell *cell, psFits *fits, pmConfig *config) 1124 1124 { 1125 1125 PS_ASSERT_PTR_NON_NULL(cell, false); 1126 1126 PS_ASSERT_FITS_NON_NULL(fits, false); 1127 1127 1128 return cellRead(cell, fits, config, FPA_READ_TYPE_ WEIGHT);1129 } 1130 1131 bool pmChipRead Weight(pmChip *chip, psFits *fits, pmConfig *config)1128 return cellRead(cell, fits, config, FPA_READ_TYPE_VARIANCE); 1129 } 1130 1131 bool pmChipReadVariance(pmChip *chip, psFits *fits, pmConfig *config) 1132 1132 { 1133 1133 PS_ASSERT_PTR_NON_NULL(chip, false); 1134 1134 PS_ASSERT_FITS_NON_NULL(fits, false); 1135 1135 1136 return chipRead(chip, fits, config, FPA_READ_TYPE_ WEIGHT);1137 } 1138 1139 bool pmFPARead Weight(pmFPA *fpa, psFits *fits, pmConfig *config)1136 return chipRead(chip, fits, config, FPA_READ_TYPE_VARIANCE); 1137 } 1138 1139 bool pmFPAReadVariance(pmFPA *fpa, psFits *fits, pmConfig *config) 1140 1140 { 1141 1141 PS_ASSERT_PTR_NON_NULL(fpa, false); 1142 1142 PS_ASSERT_FITS_NON_NULL(fits, false); 1143 1143 1144 return fpaRead(fpa, fits, config, FPA_READ_TYPE_ WEIGHT);1144 return fpaRead(fpa, fits, config, FPA_READ_TYPE_VARIANCE); 1145 1145 } 1146 1146 -
branches/pap_branch_20090128/psModules/src/camera/pmFPARead.h
r18163 r21211 4 4 * @author Paul Price, IfA 5 5 * 6 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $7 * @date $Date: 200 8-06-17 22:16:38$6 * @version $Revision: 1.15.44.1 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2009-01-29 00:33:51 $ 8 8 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii 9 9 */ … … 60 60 int numRows, // The number of rows to read 61 61 pmConfig *config 62 );62 ); 63 63 64 64 /// Return the number of readouts within a cell 65 65 /// 66 /// This function is type-independent (doesn't matter if you are interested in the image/mask/ weight).66 /// This function is type-independent (doesn't matter if you are interested in the image/mask/variance). 67 67 int pmCellNumReadouts(pmCell *cell, psFits *fits, pmConfig *config); 68 68 … … 74 74 psFits *fits, // FITS file from which to read 75 75 pmConfig *config // Configuration 76 );76 ); 77 77 78 78 /// Read a chip … … 82 82 psFits *fits, // FITS file from which to read 83 83 pmConfig *config // Configuration 84 );84 ); 85 85 86 86 /// Read an FPA … … 90 90 psFits *fits, // FITS file from which to read 91 91 pmConfig *config // Configuration 92 );92 ); 93 93 94 94 // Mask functions follow … … 126 126 psFits *fits, // FITS file from which to read 127 127 pmConfig *config // Configuration 128 );128 ); 129 129 130 130 /// Read an entire chip into the mask … … 134 134 psFits *fits, // FITS file from which to read 135 135 pmConfig *config // Configuration 136 );136 ); 137 137 138 138 /// Read an entire FPA into the mask … … 142 142 psFits *fits, // FITS file from which to read 143 143 pmConfig *config // Configuration 144 );145 146 // Weightfunctions follow147 148 /// Check to see if there is more to read when reading a readout incrementally into the weight149 bool pmReadoutMore Weight(pmReadout *readout, ///< Readout of interest150 psFits *fits,///< FITS file from which to read151 int z,///< Readout number/plane; zero-offset indexing152 int numScans,///< Number of scans (rows/cols) to read153 pmConfig *config ///< Configuration154 ); 155 156 /// Read a chunk of a readout into the weight144 ); 145 146 // Variance functions follow 147 148 /// Check to see if there is more to read when reading a readout incrementally into the variance 149 bool pmReadoutMoreVariance(pmReadout *readout, ///< Readout of interest 150 psFits *fits, ///< FITS file from which to read 151 int z, ///< Readout number/plane; zero-offset indexing 152 int numScans, ///< Number of scans (rows/cols) to read 153 pmConfig *config ///< Configuration 154 ); 155 156 /// Read a chunk of a readout into the variance 157 157 /// 158 158 /// Allows reading the readout incrementally 159 bool pmReadoutReadChunk Weight(pmReadout *readout, ///< Readout of interest160 psFits *fits, ///< FITS file from which to read161 int z, ///< Readout number/plane; zero-offset indexing162 int numScans, ///< Number of scans (rows/cols) to read163 int overlap, ///< Overlap between consecutive reads164 pmConfig *config ///< Configuration165 ); 166 167 /// Read the entire readout into the weight168 bool pmReadoutRead Weight(pmReadout *readout, ///< Readout of interest169 psFits *fits, ///< FITS file from which to read170 int z, ///< Readout number/plane; zero-offset indexing171 pmConfig *config ///< Configuration172 ); 173 174 /// Read an entire cell into the weight175 /// 176 /// Same as pmCellRead, but reads into the weightelement of the readouts.177 bool pmCellRead Weight(pmCell *cell, // Cell to read into178 psFits *fits, // FITS file from which to read179 pmConfig *config // Configuration180 );181 182 /// Read an entire chip into the weight183 /// 184 /// Same as pmChipRead, but reads into the weightelement of the readouts.185 bool pmChipRead Weight(pmChip *chip, // Chip to read into186 psFits *fits, // FITS file from which to read187 pmConfig *config // Configuration188 );189 190 /// Read an entire FPA into the weight191 /// 192 /// Same as pmFPARead, but reads into the weightelement of the readouts.193 bool pmFPARead Weight(pmFPA *fpa, // FPA to read into194 psFits *fits, // FITS file from which to read195 pmConfig *config // Configuration196 );159 bool pmReadoutReadChunkVariance(pmReadout *readout, ///< Readout of interest 160 psFits *fits, ///< FITS file from which to read 161 int z, ///< Readout number/plane; zero-offset indexing 162 int numScans, ///< Number of scans (rows/cols) to read 163 int overlap, ///< Overlap between consecutive reads 164 pmConfig *config ///< Configuration 165 ); 166 167 /// Read the entire readout into the variance 168 bool pmReadoutReadVariance(pmReadout *readout, ///< Readout of interest 169 psFits *fits, ///< FITS file from which to read 170 int z, ///< Readout number/plane; zero-offset indexing 171 pmConfig *config ///< Configuration 172 ); 173 174 /// Read an entire cell into the variance 175 /// 176 /// Same as pmCellRead, but reads into the variance element of the readouts. 177 bool pmCellReadVariance(pmCell *cell, // Cell to read into 178 psFits *fits, // FITS file from which to read 179 pmConfig *config // Configuration 180 ); 181 182 /// Read an entire chip into the variance 183 /// 184 /// Same as pmChipRead, but reads into the variance element of the readouts. 185 bool pmChipReadVariance(pmChip *chip, // Chip to read into 186 psFits *fits, // FITS file from which to read 187 pmConfig *config // Configuration 188 ); 189 190 /// Read an entire FPA into the variance 191 /// 192 /// Same as pmFPARead, but reads into the variance element of the readouts. 193 bool pmFPAReadVariance(pmFPA *fpa, // FPA to read into 194 psFits *fits, // FITS file from which to read 195 pmConfig *config // Configuration 196 ); 197 197 198 198 /// Read cell headers -
branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c
r19385 r21211 29 29 FPA_WRITE_TYPE_IMAGE, // Write image 30 30 FPA_WRITE_TYPE_MASK, // Write mask 31 FPA_WRITE_TYPE_ WEIGHT // Write weightmap31 FPA_WRITE_TYPE_VARIANCE // Write variance map 32 32 } fpaWriteType; 33 33 … … 46 46 case FPA_WRITE_TYPE_MASK: 47 47 return &hdu->masks; 48 case FPA_WRITE_TYPE_ WEIGHT:49 return &hdu-> weights;48 case FPA_WRITE_TYPE_VARIANCE: 49 return &hdu->variances; 50 50 default: 51 51 psAbort("Unknown write type: %x\n", type); … … 66 66 case FPA_WRITE_TYPE_MASK: 67 67 return pmHDUWriteMask(hdu, fits, config); 68 case FPA_WRITE_TYPE_ WEIGHT:69 return pmHDUWrite Weight(hdu, fits, config);68 case FPA_WRITE_TYPE_VARIANCE: 69 return pmHDUWriteVariance(hdu, fits, config); 70 70 default: 71 71 psAbort("Unknown write type: %x\n", type); … … 74 74 } 75 75 76 // Write a cell image/mask/ weight76 // Write a cell image/mask/variance 77 77 static bool cellWrite(pmCell *cell, // Cell to write 78 78 psFits *fits, // FITS file to which to write … … 124 124 } 125 125 126 // Write a chip image/mask/ weight126 // Write a chip image/mask/variance 127 127 static bool chipWrite(pmChip *chip, // Chip to write 128 128 psFits *fits, // FITS file to which to write … … 188 188 189 189 190 // Write an FPA image/mask/ weight190 // Write an FPA image/mask/variance 191 191 static bool fpaWrite(pmFPA *fpa, // FPA to write 192 192 psFits *fits, // FITS file to which to write … … 417 417 418 418 419 bool pmCellWrite Weight(pmCell *cell, psFits *fits, pmConfig *config, bool blank)419 bool pmCellWriteVariance(pmCell *cell, psFits *fits, pmConfig *config, bool blank) 420 420 { 421 421 PS_ASSERT_PTR_NON_NULL(cell, false); 422 422 PS_ASSERT_PTR_NON_NULL(fits, false); 423 return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_ WEIGHT);424 } 425 426 bool pmChipWrite Weight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)423 return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_VARIANCE); 424 } 425 426 bool pmChipWriteVariance(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) 427 427 { 428 428 PS_ASSERT_PTR_NON_NULL(chip, false); 429 429 PS_ASSERT_PTR_NON_NULL(fits, false); 430 return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_ WEIGHT);431 } 432 433 bool pmFPAWrite Weight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)430 return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE); 431 } 432 433 bool pmFPAWriteVariance(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) 434 434 { 435 435 PS_ASSERT_PTR_NON_NULL(fpa, false); 436 436 PS_ASSERT_PTR_NON_NULL(fits, false); 437 return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_ WEIGHT);437 return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE); 438 438 } 439 439 -
branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.h
r18163 r21211 4 4 * @author Paul Price, IfA 5 5 * 6 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $7 * @date $Date: 200 8-06-17 22:16:38$6 * @version $Revision: 1.11.44.1 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2009-01-29 00:33:51 $ 8 8 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii 9 9 */ … … 24 24 psFits *fits, ///< FITS file to which to write 25 25 int z ///< Image plane to write 26 );26 ); 27 27 28 28 /// Write a cell to a FITS file … … 36 36 pmConfig *config, ///< Configuration 37 37 bool blank ///< Write a blank PHU? 38 );38 ); 39 39 40 40 /// Write a chip to a FITS file … … 49 49 bool blank, ///< Write a blank PHU? 50 50 bool recurse ///< Recurse to lower levels? 51 );51 ); 52 52 53 53 /// Write an FPA to a FITS file … … 62 62 bool blank, ///< Write a blank PHU? 63 63 bool recurse ///< Recurse to lower levels? 64 );64 ); 65 65 66 66 /// Write a cell mask to a FITS file … … 74 74 pmConfig *config, ///< Configuration 75 75 bool blank ///< Write a blank PHU? 76 );76 ); 77 77 78 78 /// Write a chip mask to a FITS file … … 88 88 bool blank, ///< Write a blank PHU? 89 89 bool recurse ///< Recurse to lower levels? 90 );90 ); 91 91 92 92 /// Write an FPA mask to a FITS file … … 102 102 bool blank, ///< Write a blank PHU? 103 103 bool recurse ///< Recurse to lower levels? 104 );104 ); 105 105 106 /// Write a cell weightto a FITS file106 /// Write a cell variance to a FITS file 107 107 /// 108 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weightpixels if required. A blank (i.e., image-less108 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less 109 109 /// header) is written only if specifically requested. Writes the concepts to the various locations, and then 110 /// the HDU weight to the FITS file. This function should be called at the beginning of the output cell loop111 /// with blank=true in order to produce the correct file structure.112 bool pmCellWrite Weight(pmCell *cell, ///< Cell to write113 psFits *fits, ///< FITS file to which to write114 pmConfig *config, ///< Configuration115 bool blank ///< Write a blank PHU?116 );110 /// the HDU variance to the FITS file. This function should be called at the beginning of the output cell 111 /// loop with blank=true in order to produce the correct file structure. 112 bool pmCellWriteVariance(pmCell *cell, ///< Cell to write 113 psFits *fits, ///< FITS file to which to write 114 pmConfig *config, ///< Configuration 115 bool blank ///< Write a blank PHU? 116 ); 117 117 118 /// Write a chip weightto a FITS file118 /// Write a chip variance to a FITS file 119 119 /// 120 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weightpixels if required. A blank (i.e., image-less120 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less 121 121 /// header) is written only if specifically requested. Writes the concepts to the various locations, and then 122 /// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at123 /// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct122 /// the HDU variance to the FITS file, optionally recursing to lower levels. This function should be called 123 /// at the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct 124 124 /// file structure. 125 bool pmChipWrite Weight(pmChip *chip, ///< Chip to write126 psFits *fits, ///< FITS file to which to write127 pmConfig *config, ///< Configuration128 bool blank, ///< Write a blank PHU?129 bool recurse ///< Recurse to lower levels?130 );125 bool pmChipWriteVariance(pmChip *chip, ///< Chip to write 126 psFits *fits, ///< FITS file to which to write 127 pmConfig *config, ///< Configuration 128 bool blank, ///< Write a blank PHU? 129 bool recurse ///< Recurse to lower levels? 130 ); 131 131 132 /// Write an FPA weightto a FITS file132 /// Write an FPA variance to a FITS file 133 133 /// 134 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weightpixels if required. A blank (i.e., image-less134 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required. A blank (i.e., image-less 135 135 /// header) is written only if specifically requested. Writes the concepts to the various locations, and then 136 /// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at137 /// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct136 /// the HDU variance to the FITS file, optionally recursing to lower levels. This function should be called 137 /// at the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct 138 138 /// file structure. 139 bool pmFPAWrite Weight(pmFPA *fpa, ///< FPA to write140 psFits *fits, ///< FITS file to which to write141 pmConfig *config, ///< Configuration142 bool blank, ///< Write a blank PHU?143 bool recurse ///< Recurse to lower levels?144 );139 bool pmFPAWriteVariance(pmFPA *fpa, ///< FPA to write 140 psFits *fits, ///< FITS file to which to write 141 pmConfig *config, ///< Configuration 142 bool blank, ///< Write a blank PHU? 143 bool recurse ///< Recurse to lower levels? 144 ); 145 145 146 146 … … 153 153 const pmCell *cell, ///< Cell containing FITS table in the analysis metadata 154 154 const char *name ///< Name for the table data, and the extension name 155 );155 ); 156 156 157 157 int pmChipWriteTable(psFits *fits, ///< FITS file to which to write 158 158 const pmChip *chip, ///< Chip containing cells with tables to write 159 159 const char *name ///< Name for the table data, and the extension name 160 );160 ); 161 161 162 162 int pmFPAWriteTable(psFits *fits, ///< FITS file to which to write 163 163 const pmFPA *fpa, ///< FPA containing cells with tables to write 164 164 const char *name ///< Name for the table data, and the extension name 165 );165 ); 166 166 167 167 // XXX better name, please -
branches/pap_branch_20090128/psModules/src/camera/pmFPAfile.c
r19307 r21211 473 473 return PM_FPA_FILE_MASK; 474 474 } 475 if (!strcasecmp(type, " WEIGHT")) {476 return PM_FPA_FILE_ WEIGHT;475 if (!strcasecmp(type, "VARIANCE")) { 476 return PM_FPA_FILE_VARIANCE; 477 477 } 478 478 if (!strcasecmp(type, "FRINGE")) { … … 527 527 case PM_FPA_FILE_MASK: 528 528 return ("MASK"); 529 case PM_FPA_FILE_ WEIGHT:530 return (" WEIGHT");529 case PM_FPA_FILE_VARIANCE: 530 return ("VARIANCE"); 531 531 case PM_FPA_FILE_FRINGE: 532 532 return ("FRINGE"); -
branches/pap_branch_20090128/psModules/src/camera/pmFPAfile.h
r19307 r21211 4 4 * @author EAM, IfA 5 5 * 6 * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $7 * @date $Date: 200 8-09-02 19:06:16$6 * @version $Revision: 1.33.22.1 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2009-01-29 00:33:51 $ 8 8 * Copyright 2004-2005 Institute for Astronomy, University of Hawaii 9 9 */ … … 37 37 PM_FPA_FILE_IMAGE, 38 38 PM_FPA_FILE_MASK, 39 PM_FPA_FILE_ WEIGHT,39 PM_FPA_FILE_VARIANCE, 40 40 PM_FPA_FILE_FRINGE, 41 41 PM_FPA_FILE_DARK, -
branches/pap_branch_20090128/psModules/src/camera/pmFPAfileFitsIO.c
r18601 r21211 148 148 case PM_FPA_FILE_IMAGE: 149 149 case PM_FPA_FILE_MASK: 150 case PM_FPA_FILE_ WEIGHT:150 case PM_FPA_FILE_VARIANCE: 151 151 case PM_FPA_FILE_HEADER: 152 152 case PM_FPA_FILE_FRINGE: … … 329 329 } 330 330 331 bool pmFPAviewReadFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) 332 { 333 PS_ASSERT_PTR_NON_NULL(view, false); 334 PS_ASSERT_PTR_NON_NULL(file, false); 335 return fpaViewReadFitsImage(view, file, config, pmFPAReadWeight, pmChipReadWeight, pmCellReadWeight); 331 bool pmFPAviewReadFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config) 332 { 333 PS_ASSERT_PTR_NON_NULL(view, false); 334 PS_ASSERT_PTR_NON_NULL(file, false); 335 return fpaViewReadFitsImage(view, file, config, pmFPAReadVariance, pmChipReadVariance, 336 pmCellReadVariance); 336 337 } 337 338 … … 347 348 PS_ASSERT_PTR_NON_NULL(view, false); 348 349 PS_ASSERT_PTR_NON_NULL(file, false); 349 return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, pmCellReadHeaderSet); 350 return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, 351 pmCellReadHeaderSet); 350 352 } 351 353 … … 432 434 } 433 435 434 bool pmFPAviewWriteFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) 435 { 436 PS_ASSERT_PTR_NON_NULL(view, false); 437 PS_ASSERT_PTR_NON_NULL(file, false); 438 return fpaViewWriteFitsImage(view, file, config, pmFPAWriteWeight, pmChipWriteWeight, pmCellWriteWeight); 436 bool pmFPAviewWriteFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config) 437 { 438 PS_ASSERT_PTR_NON_NULL(view, false); 439 PS_ASSERT_PTR_NON_NULL(file, false); 440 return fpaViewWriteFitsImage(view, file, config, pmFPAWriteVariance, pmChipWriteVariance, 441 pmCellWriteVariance); 439 442 } 440 443 -
branches/pap_branch_20090128/psModules/src/camera/pmFPAfileFitsIO.h
r18601 r21211 5 5 * @author PAP, IfA 6 6 * 7 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $8 * @date $Date: 200 8-07-17 22:38:15$7 * @version $Revision: 1.15.36.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2009-01-29 00:33:51 $ 9 9 * Copyright 2004-2005 Institute for Astronomy, University of Hawaii 10 10 */ … … 27 27 pmConfig *config 28 28 ); 29 /// Read a weightmap into the current view30 bool pmFPAviewReadFits Weight(const pmFPAview *view, ///< View specifying level of interest29 /// Read a variance map into the current view 30 bool pmFPAviewReadFitsVariance(const pmFPAview *view, ///< View specifying level of interest 31 31 pmFPAfile *file, ///< FPA file into which to read 32 32 pmConfig *config … … 57 57 ); 58 58 59 /// Write the weightmap for the specified view60 bool pmFPAviewWriteFits Weight(const pmFPAview *view, ///< View specifying level of interest59 /// Write the variance map for the specified view 60 bool pmFPAviewWriteFitsVariance(const pmFPAview *view, ///< View specifying level of interest 61 61 pmFPAfile *file, ///< FPA file to write 62 62 pmConfig *config ///< Configuration -
branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c
r20937 r21211 173 173 status = pmFPAviewReadFitsMask(view, file, config); 174 174 break; 175 case PM_FPA_FILE_ WEIGHT:176 status = pmFPAviewReadFits Weight(view, file, config);175 case PM_FPA_FILE_VARIANCE: 176 status = pmFPAviewReadFitsVariance(view, file, config); 177 177 break; 178 178 case PM_FPA_FILE_HEADER: … … 262 262 case PM_FPA_FILE_IMAGE: 263 263 case PM_FPA_FILE_MASK: 264 case PM_FPA_FILE_ WEIGHT:264 case PM_FPA_FILE_VARIANCE: 265 265 case PM_FPA_FILE_FRINGE: 266 266 case PM_FPA_FILE_DARK: { … … 427 427 status = pmFPAviewWriteFitsMask(view, file, config); 428 428 break; 429 case PM_FPA_FILE_ WEIGHT:430 status = pmFPAviewWriteFits Weight(view, file, config);429 case PM_FPA_FILE_VARIANCE: 430 status = pmFPAviewWriteFitsVariance(view, file, config); 431 431 break; 432 432 case PM_FPA_FILE_HEADER: … … 524 524 case PM_FPA_FILE_IMAGE: 525 525 case PM_FPA_FILE_MASK: 526 case PM_FPA_FILE_ WEIGHT:526 case PM_FPA_FILE_VARIANCE: 527 527 case PM_FPA_FILE_HEADER: 528 528 case PM_FPA_FILE_FRINGE: … … 589 589 case PM_FPA_FILE_IMAGE: 590 590 case PM_FPA_FILE_MASK: 591 case PM_FPA_FILE_ WEIGHT:591 case PM_FPA_FILE_VARIANCE: 592 592 case PM_FPA_FILE_HEADER: 593 593 case PM_FPA_FILE_FRINGE: … … 730 730 psString tmpName = pmConfigConvertFilename (file->filename, config, create, false); 731 731 if (!tmpName) { 732 psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename);733 return false;732 psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename); 733 return false; 734 734 } 735 735 psFree (file->filename); … … 740 740 case PM_FPA_FILE_IMAGE: 741 741 case PM_FPA_FILE_MASK: 742 case PM_FPA_FILE_ WEIGHT:742 case PM_FPA_FILE_VARIANCE: 743 743 case PM_FPA_FILE_HEADER: 744 744 case PM_FPA_FILE_FRINGE: … … 783 783 // determine camera if not specified already 784 784 // XXX can I actually reach this with camera not specified?? 785 psMetadata *camera = NULL;786 psString formatName = NULL;787 psString cameraName = NULL;785 psMetadata *camera = NULL; 786 psString formatName = NULL; 787 psString cameraName = NULL; 788 788 file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true); 789 789 if (!file->format) { … … 794 794 psFree(phu); 795 795 796 pmFPA *newFPA = pmFPAConstruct (camera, formatName);797 if (!newFPA) {798 psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName);799 psFree(camera);800 psFree(formatName);801 return NULL;802 }803 psFree(camera);804 psFree(formatName);805 psFree(cameraName);806 807 // XXX this is really dangerous...808 psFree (file->fpa);809 file->fpa = newFPA;796 pmFPA *newFPA = pmFPAConstruct (camera, formatName); 797 if (!newFPA) { 798 psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName); 799 psFree(camera); 800 psFree(formatName); 801 return NULL; 802 } 803 psFree(camera); 804 psFree(formatName); 805 psFree(cameraName); 806 807 // XXX this is really dangerous... 808 psFree (file->fpa); 809 file->fpa = newFPA; 810 810 } 811 811 … … 938 938 case PM_FPA_FILE_IMAGE: 939 939 case PM_FPA_FILE_MASK: 940 case PM_FPA_FILE_ WEIGHT:940 case PM_FPA_FILE_VARIANCE: 941 941 case PM_FPA_FILE_DARK: 942 942 case PM_FPA_FILE_FRINGE: -
branches/pap_branch_20090128/psModules/src/camera/pmHDU.c
r21183 r21211 45 45 psFree(hdu->header); 46 46 psFree(hdu->images); 47 psFree(hdu-> weights);47 psFree(hdu->variances); 48 48 psFree(hdu->masks); 49 49 } … … 69 69 hdu->header = NULL; 70 70 hdu->images = NULL; 71 hdu-> weights = NULL;71 hdu->variances = NULL; 72 72 hdu->masks = NULL; 73 73 … … 155 155 } 156 156 157 bool pmHDURead Weight(pmHDU *hdu, psFits *fits)158 { 159 PS_ASSERT_PTR_NON_NULL(hdu, false); 160 PS_ASSERT_PTR_NON_NULL(fits, false); 161 162 return hduRead(hdu, &hdu-> weights, fits);157 bool pmHDUReadVariance(pmHDU *hdu, psFits *fits) 158 { 159 PS_ASSERT_PTR_NON_NULL(hdu, false); 160 PS_ASSERT_PTR_NON_NULL(fits, false); 161 162 return hduRead(hdu, &hdu->variances, fits); 163 163 } 164 164 … … 235 235 } 236 236 237 bool pmHDUWrite Weight(pmHDU *hdu, psFits *fits, const pmConfig *config)237 bool pmHDUWriteVariance(pmHDU *hdu, psFits *fits, const pmConfig *config) 238 238 { 239 239 PS_ASSERT_PTR_NON_NULL(hdu, false); … … 241 241 242 242 psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask 243 return hduWrite(hdu, hdu-> weights, hdu->masks, maskVal, fits);244 } 243 return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits); 244 } -
branches/pap_branch_20090128/psModules/src/camera/pmHDU.h
r19385 r21211 4 4 * @author Paul Price, IfA 5 5 * 6 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $7 * @date $Date: 200 8-09-05 08:21:35$6 * @version $Revision: 1.9.22.1 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2009-01-29 00:33:51 $ 8 8 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii 9 9 */ … … 20 20 /// An instance of the FITS Header Data Unit 21 21 /// 22 /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weightdata, but these are22 /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and variance data, but these are 23 23 /// stored here for convenience --- it keeps all the relevant data about the image in one place. 24 24 typedef struct … … 29 29 psMetadata *header; ///< The FITS header, or NULL if primary for FITS; or section info 30 30 psArray *images; ///< The pixel data 31 psArray * weights; ///< The pixel data31 psArray *variances; ///< The pixel data 32 32 psArray *masks; ///< The pixel data 33 33 } … … 60 60 ); 61 61 62 /// Read the HDU header and weightmap62 /// Read the HDU header and variance map 63 63 /// 64 64 /// Moves to the appropriate extension 65 bool pmHDURead Weight(pmHDU *hdu, ///< HDU to read66 psFits *fits ///< FITS file to read from65 bool pmHDUReadVariance(pmHDU *hdu, ///< HDU to read 66 psFits *fits ///< FITS file to read from 67 67 ); 68 68 … … 79 79 ); 80 80 81 /// Write the HDU header and weightmap82 bool pmHDUWrite Weight(pmHDU *hdu, ///< HDU to write83 psFits *fits, ///< FITS file to write to84 const pmConfig *config ///< Configuration81 /// Write the HDU header and variance map 82 bool pmHDUWriteVariance(pmHDU *hdu, ///< HDU to write 83 psFits *fits, ///< FITS file to write to 84 const pmConfig *config ///< Configuration 85 85 ); 86 86 -
branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c
r18227 r21211 274 274 pmReadout *readout = cell->readouts->data[0]; // The first readout, as representative 275 275 // The proper image, used to get the size 276 psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout-> weight);276 psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout->variance); 277 277 if (!image) { 278 278 continue; … … 283 283 image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows); 284 284 } 285 if (readout-> weight&&286 (readout-> weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) {287 psWarning("Image and weighthave different sizes (%dx%d vs %dx%d)!\n",288 image->numCols, image->numRows, readout-> weight->numCols, readout->weight->numRows);285 if (readout->variance && 286 (readout->variance->numCols != image->numCols || readout->variance->numRows != image->numRows)) { 287 psWarning("Image and variance have different sizes (%dx%d vs %dx%d)!\n", 288 image->numCols, image->numRows, readout->variance->numCols, readout->variance->numRows); 289 289 } 290 290 // New reference … … 355 355 psElemType imageType = 0; // Type of readout images 356 356 psElemType maskType = 0; // Type of readout masks 357 psElemType weightType = 0; // Type of readout weights357 psElemType varianceType = 0; // Type of readout variances 358 358 { 359 359 psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells … … 381 381 maskType = checkTypes(maskType, readout->mask->type.type); 382 382 } 383 if (!hdu-> weights && readout->weight) {384 weightType = checkTypes(weightType, readout->weight->type.type);383 if (!hdu->variances && readout->variance) { 384 varianceType = checkTypes(varianceType, readout->variance->type.type); 385 385 } 386 386 } … … 388 388 psFree(iter); 389 389 } 390 if (numReadouts == 0 || (imageType == 0 && maskType == 0 && weightType == 0)) {390 if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) { 391 391 // Nothing from which to create an HDU 392 392 psFree(cells); … … 421 421 } 422 422 } 423 if ( weightType) {424 hdu-> weights = psArrayAlloc(numReadouts);423 if (varianceType) { 424 hdu->variances = psArrayAlloc(numReadouts); 425 425 for (int i = 0; i < numReadouts; i++) { 426 psImage * weight = psImageAlloc(xSize, ySize, weightType);427 psImageInit( weight, 0.0);428 hdu-> weights->data[i] = weight;426 psImage *variance = psImageAlloc(xSize, ySize, varianceType); 427 psImageInit(variance, 0.0); 428 hdu->variances->data[i] = variance; 429 429 } 430 430 } … … 450 450 psArray *hduImages = hdu->images; // Array of images in the HDU 451 451 psArray *hduMasks = hdu->masks; // Array of masks in the HDU 452 psArray *hdu Weights = hdu->weights; // Array of weights in the HDU452 psArray *hduVariances = hdu->variances; // Array of variances in the HDU 453 453 for (int i = 0; i < readouts->n; i++) { 454 454 pmReadout *readout = readouts->data[i]; // The readout of interest … … 467 467 readout->mask = new; 468 468 } 469 if (readout-> weight) {470 psImage *new = pasteImage(hdu Weights->data[i], readout->weight, trimsec);471 psFree(readout-> weight);472 readout-> weight= new;469 if (readout->variance) { 470 psImage *new = pasteImage(hduVariances->data[i], readout->variance, trimsec); 471 psFree(readout->variance); 472 readout->variance = new; 473 473 } 474 474 … … 581 581 return false; 582 582 } 583 if (hdu->images && hdu->masks && hdu-> weights) {583 if (hdu->images && hdu->masks && hdu->variances) { 584 584 // It's already here! 585 585 return true; … … 631 631 return generateForCells(chip); 632 632 } 633 if (hdu->images && hdu->masks && hdu-> weights) {633 if (hdu->images && hdu->masks && hdu->variances) { 634 634 // It's already here! 635 635 return true; … … 680 680 return generateForChips(fpa); 681 681 } 682 if (hdu->images && hdu->masks && hdu-> weights) {682 if (hdu->images && hdu->masks && hdu->variances) { 683 683 // It's already here! 684 684 return true; -
branches/pap_branch_20090128/psModules/src/camera/pmHDUUtils.c
r18554 r21211 18 18 19 19 for (int i = 0; i < fpa->chips->n; i++) { 20 pmChip *chip = fpa->chips->data[i];21 if (!chip) continue;22 if (chip->hdu) return chip->hdu;23 if (!chip->cells) continue;24 for (int j = 0; j < chip->cells->n; j++) {25 pmCell *cell = chip->cells->data[j];26 if (!cell) continue;27 if (cell->hdu) return cell->hdu;28 }20 pmChip *chip = fpa->chips->data[i]; 21 if (!chip) continue; 22 if (chip->hdu) return chip->hdu; 23 if (!chip->cells) continue; 24 for (int j = 0; j < chip->cells->n; j++) { 25 pmCell *cell = chip->cells->data[j]; 26 if (!cell) continue; 27 if (cell->hdu) return cell->hdu; 28 } 29 29 } 30 30 return NULL; … … 158 158 159 159 INDENT(fd, level + 1); 160 if (hdu-> weights) {161 fprintf(fd, " Weights:\n");162 for (long i = 0; i < hdu-> weights->n; i++) {163 psImage * weight = hdu->weights->data[i]; // Weightimage of interest160 if (hdu->variances) { 161 fprintf(fd, "Variances:\n"); 162 for (long i = 0; i < hdu->variances->n; i++) { 163 psImage *variance = hdu->variances->data[i]; // Variance image of interest 164 164 INDENT(fd, level + 2); 165 fprintf(fd, "%ld: %dx%d\n", i, weight->numCols, weight->numRows);165 fprintf(fd, "%ld: %dx%d\n", i, variance->numCols, variance->numRows); 166 166 } 167 167 } else { 168 fprintf(fd, "NO weights.\n");168 fprintf(fd, "NO variances.\n"); 169 169 } 170 170 -
branches/pap_branch_20090128/psModules/src/camera/pmReadoutStack.c
r21183 r21211 12 12 // XXX should it be an error for the image to exist? 13 13 psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image 14 const char *name, // Name of image in analysis metadata15 int numCols, int numRows, // Expected size of image16 psElemType type, // Expected type of image17 double init // Initial value14 const char *name, // Name of image in analysis metadata 15 int numCols, int numRows, // Expected size of image 16 psElemType type, // Expected type of image 17 double init // Initial value 18 18 ) 19 19 { … … 24 24 25 25 if (!psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image)) { 26 psAbort ("analysis image already exists");26 psAbort ("analysis image already exists"); 27 27 } 28 28 psImageInit(image, init); … … 35 35 // XXX not sure why this should call psMemIncrRefCounter 36 36 psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image 37 const char *name // Name of image in analysis metadata37 const char *name // Name of image in analysis metadata 38 38 ) 39 39 { … … 78 78 79 79 // XXX for the moment, use col0, row0, numCols, numRows supplied from the outside 80 bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool weight, psImageMaskType blank)80 bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool variance, psImageMaskType blank) 81 81 { 82 82 PS_ASSERT_PTR_NON_NULL(readout, false); … … 92 92 93 93 if (mask) { 94 // XXX is this an error?94 // XXX is this an error? 95 95 if (readout->mask) return false; 96 readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);97 psImageInit(readout->mask, blank);98 } 99 100 if ( weight) {101 // XXX is this an error?102 if (readout-> weight) return false;103 readout->weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);104 psImageInit(readout->weight, NAN);96 readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); 97 psImageInit(readout->mask, blank); 98 } 99 100 if (variance) { 101 // XXX is this an error? 102 if (readout->variance) return false; 103 readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 104 psImageInit(readout->variance, NAN); 105 105 } 106 106 … … 162 162 163 163 bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, 164 int numCols, int numRows, bool mask, bool weight,164 int numCols, int numRows, bool mask, bool variance, 165 165 psImageMaskType blank) 166 166 { … … 203 203 } 204 204 205 if ( weight) {206 if (!readout-> weight) {207 readout-> weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);208 psImageInit(readout-> weight, NAN);209 } 210 if (readout-> weight->numCols < numCols || readout->weight->numRows < numRows) {211 psImage *new Weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);212 psImageInit(new Weight, NAN);213 psImageOverlaySection(new Weight, readout->weight, readout->col0, readout->row0, "=");214 psFree(readout-> weight);215 readout-> weight = newWeight;205 if (variance) { 206 if (!readout->variance) { 207 readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 208 psImageInit(readout->variance, NAN); 209 } 210 if (readout->variance->numCols < numCols || readout->variance->numRows < numRows) { 211 psImage *newVariance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 212 psImageInit(newVariance, NAN); 213 psImageOverlaySection(newVariance, readout->variance, readout->col0, readout->row0, "="); 214 psFree(readout->variance); 215 readout->variance = newVariance; 216 216 } 217 217 } -
branches/pap_branch_20090128/psModules/src/detrend/pmShutterCorrection.c
r21183 r21211 383 383 psArray *images = psArrayAlloc(num);// Array of images 384 384 psArray *masks = NULL; // Array of masks 385 psArray * weights = NULL; // Array of weights385 psArray *variances = NULL; // Array of variances 386 386 psVector *exptimes = psVectorAlloc(num, PS_TYPE_F32); // Vector of exposure times 387 387 … … 392 392 masks = psArrayAlloc(num); 393 393 } 394 if (readout-> weight)394 if (readout->variance) 395 395 { 396 weights = psArrayAlloc(num);396 variances = psArrayAlloc(num); 397 397 } 398 398 } … … 472 472 } 473 473 474 psImage * weight = readout->weight; // Weightmap of interest475 if ( weight) {476 if (! weights) {477 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have weights.\n");474 psImage *variance = readout->variance; // Variance map of interest 475 if (variance) { 476 if (!variances) { 477 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have variances.\n"); 478 478 goto MEASURE_ERROR; 479 479 } 480 weights->data[i] = psMemIncrRefCounter(weight);481 482 if ( weight->type.type != PS_TYPE_F32) {483 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for weights: %x\n", weight->type.type);480 variances->data[i] = psMemIncrRefCounter(variance); 481 482 if (variance->type.type != PS_TYPE_F32) { 483 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for variances: %x\n", variance->type.type); 484 484 goto MEASURE_ERROR; 485 485 } 486 if ( weight->numRows != numRows || weight->numCols != numCols) {486 if (variance->numRows != numRows || variance->numCols != numCols) { 487 487 psError(PS_ERR_BAD_PARAMETER_SIZE, true, 488 " Weight sizes don't match: %dx%d vs %dx%d\n", weight->numCols, weight->numRows,488 "Variance sizes don't match: %dx%d vs %dx%d\n", variance->numCols, variance->numRows, 489 489 numCols, numRows); 490 490 goto MEASURE_ERROR; 491 491 } 492 } else if ( weights) {493 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have weights.\n");492 } else if (variances) { 493 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not all readouts have variances.\n"); 494 494 goto MEASURE_ERROR; 495 495 } … … 566 566 psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef); 567 567 568 // Check the weights569 if ( weights && nIter > 1) {570 for (int i = 0; i < weights->n && nIter > 1; i++) {571 psImage * weight = weights->data[i]; // Weightimage572 if (! weight) {573 // We don't have weights, so no realistic errors: turn off iteration568 // Check the variances 569 if (variances && nIter > 1) { 570 for (int i = 0; i < variances->n && nIter > 1; i++) { 571 psImage *variance = variances->data[i]; // Variance image 572 if (!variance) { 573 // We don't have variances, so no realistic errors: turn off iteration 574 574 if (nIter > 0) { 575 psWarning("Not all images have weights --- turning iteration off.\n");575 psWarning("Not all images have variances --- turning iteration off.\n"); 576 576 } 577 577 nIter = 1; … … 594 594 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal); 595 595 } 596 psImage * weight; // Weightimage597 if ( weights && (weight = weights->data[i])) {598 errors->data.F32[i] = sqrtf( weight->data.F32[y][x]) * refs->data.F32[i];596 psImage *variance; // Variance image 597 if (variances && (variance = variances->data[i])) { 598 errors->data.F32[i] = sqrtf(variance->data.F32[y][x]) * refs->data.F32[i]; 599 599 } else { 600 600 errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i]; … … 645 645 psFree(images); 646 646 psFree(masks); 647 psFree( weights);647 psFree(variances); 648 648 psFree(refs); 649 649 psFree(regions); … … 878 878 PS_ASSERT_IMAGE_SIZE(readout->mask, data->numCols, data->numRows, NULL); 879 879 } 880 if (readout-> weight) {881 PS_ASSERT_IMAGE_NON_NULL(readout-> weight, NULL);882 PS_ASSERT_IMAGE_TYPE(readout-> weight, PS_TYPE_F32, NULL);883 PS_ASSERT_IMAGE_SIZE(readout-> weight, data->numCols, data->numRows, NULL);880 if (readout->variance) { 881 PS_ASSERT_IMAGE_NON_NULL(readout->variance, NULL); 882 PS_ASSERT_IMAGE_TYPE(readout->variance, PS_TYPE_F32, NULL); 883 PS_ASSERT_IMAGE_SIZE(readout->variance, data->numCols, data->numRows, NULL); 884 884 } 885 885 … … 1133 1133 mask->data.PS_TYPE_VECTOR_MASK_DATA[r] = (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal); 1134 1134 } 1135 if (readout-> weight) {1136 errors->data.F32[r] = sqrtf(readout-> weight->data.F32[yIn][xIn]) * ref;1135 if (readout->variance) { 1136 errors->data.F32[r] = sqrtf(readout->variance->data.F32[yIn][xIn]) * ref; 1137 1137 } else { 1138 1138 // XXX guess that the input data is Poisson distributed; if we go negative, force high -
branches/pap_branch_20090128/psModules/src/imcombine/pmPSFEnvelope.c
r21183 r21211 36 36 #define PEAK_FLUX 1.0e4 // Peak flux for each source 37 37 #define SKY_VALUE 0.0e0 // Sky value for fake image 38 #define WEIGHT_VAL 3.0 // Weightingfor image39 #define WEIGHT_FACTOR 10.0 // Factor to multiply image by to get weighting38 #define VARIANCE_VAL 3.0 // Variance for image 39 #define VARIANCE_FACTOR 10.0 // Factor to multiply image by to get variance 40 40 #define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF 41 41 #define SOURCE_FIT_ITERATIONS 100 // Number of iterations for source fitting … … 149 149 pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source 150 150 psAssert (model, "failed to generate model: should this be an error or not?"); 151 float srcRadius = model->modelRadius(model->params, PS_SQR( WEIGHT_VAL)); // Radius for source151 float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source 152 152 if (srcRadius > maxRadius) { 153 153 maxRadius = srcRadius; … … 216 216 psFree(envelope); 217 217 218 // XXX Setting the weightseems to be an art218 // XXX Setting the variance seems to be an art 219 219 // Can't set it too high so that pixels are rejected as insignificant 220 220 // Can't set it too low so that it's hard to get to the minimum 221 221 // Have also tried: 222 // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); 223 // *** readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(WEIGHT_FACTOR, PS_TYPE_F32)); 224 readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(WEIGHT_VAL, PS_TYPE_F32)); 222 // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", readout->image); 223 // *** readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "*", psScalarAlloc(VARIANCE_FACTOR, PS_TYPE_F32)); 224 readout->variance = (psImage*)psBinaryOp(NULL, readout->image, "+", 225 psScalarAlloc(VARIANCE_VAL, PS_TYPE_F32)); 225 226 readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); 226 227 psImageInit(readout->mask, 0); … … 242 243 243 244 psFree(source->pixels); 244 psFree(source-> weight);245 psFree(source->variance); 245 246 psFree(source->maskView); 246 247 psFree(source->maskObj); 247 248 source->pixels = NULL; 248 source-> weight= NULL;249 source->variance = NULL; 249 250 source->maskView = NULL; 250 251 source->maskObj = NULL; … … 280 281 options->psfFieldYo = 0; 281 282 282 pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, WEIGHT_VAL, true);283 pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); 283 284 284 285 pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff); -
branches/pap_branch_20090128/psModules/src/imcombine/pmReadoutCombine.c
r21183 r21211 37 37 params->iter = 1; 38 38 params->rej = INFINITY; 39 params-> weights = false;39 params->variances = false; 40 40 41 41 return params; … … 75 75 psStringAppend(&comment, "Combining using statistic: %x", params->combine); 76 76 if (!hdu->header) { 77 hdu->header = psMetadataAlloc();77 hdu->header = psMetadataAlloc(); 78 78 } 79 79 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); … … 82 82 // note the clipping parameters, if used 83 83 if (params->combine == PS_STAT_CLIPPED_MEAN) { 84 psString comment = NULL; // Comment to add to header85 psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej);86 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");87 psFree(comment);88 } 89 90 // note the use of weights91 if (params-> weights) {84 psString comment = NULL; // Comment to add to header 85 psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); 86 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); 87 psFree(comment); 88 } 89 90 // note the use of variances 91 if (params->variances) { 92 92 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, 93 "Using input weights to combine images", "");93 "Using input variances to combine images", ""); 94 94 } 95 95 … … 121 121 122 122 // generate the required output images based on the specified sizes 123 pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params-> weights, params->blank);123 pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->variances, params->blank); 124 124 psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); 125 125 … … 177 177 return false; 178 178 } 179 if (params-> weights && !readout->weight) {179 if (params->variances && !readout->variance) { 180 180 psError(PS_ERR_UNEXPECTED_NULL, true, 181 "Rejection based on weights requested, but no weights supplied for image %d.\n", i);181 "Rejection based on variances requested, but no variances supplied for image %d.\n", i); 182 182 return false; 183 183 } … … 190 190 } 191 191 192 // pthread_t id = pthread_self(); 193 // char name[64]; 194 // sprintf (name, "%x", (unsigned int) id); 195 // psTimerStart (name); 196 197 psStatsOptions combineStdev = 0; // Statistics option for weights 192 #if 0 193 pthread_t id = pthread_self(); 194 char name[64]; 195 sprintf(name, "%x", (unsigned int)id); 196 psTimerStart(name); 197 #endif 198 199 psStatsOptions combineStdev = 0; // Statistics option for variances 198 200 switch (params->combine) { 199 201 case PS_STAT_SAMPLE_MEAN: … … 250 252 psVectorMaskType *maskData = mask->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference mask 251 253 252 psVector * weights = NULL; // Stack of weights253 psVector *errors = NULL; // Stack of errors (sqrt of variance /weights), for psVectorStats254 psF32 * weightsData = NULL; // Dereference weights255 if (params-> weights) {256 weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights257 weightsData = weights->data.F32;254 psVector *variances = NULL; // Stack of variances 255 psVector *errors = NULL; // Stack of errors (sqrt of variance), for psVectorStats 256 psF32 *variancesData = NULL; // Dereference variances 257 if (params->variances) { 258 variances = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of variances 259 variancesData = variances->data.F32; 258 260 } 259 261 psVector *index = NULL; // The indices to sort the pixels … … 279 281 psF32 **outputImage = output->image->data.F32; // Output image 280 282 psImageMaskType **outputMask = output->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Output mask 281 psF32 **output Weight = NULL; // Output weightmap282 if (output-> weight) {283 output Weight = output->weight->data.F32;283 psF32 **outputVariance = NULL; // Output variance map 284 if (output->variance) { 285 outputVariance = output->variance->data.F32; 284 286 } 285 287 … … 323 325 } 324 326 325 if (params-> weights) {326 weightsData[r] = readout->weight->data.F32[yIn][xIn];327 if (params->variances) { 328 variancesData[r] = readout->variance->data.F32[yIn][xIn]; 327 329 } 328 330 … … 332 334 if (scale) { 333 335 pixelsData[r] *= invScale->data.F32[r]; 334 if (params-> weights) {335 weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r];336 if (params->variances) { 337 variancesData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; 336 338 } 337 339 } … … 376 378 377 379 // XXXXX this step probably is very expensive : convert errors to variance everywhere? 378 if (params-> weights) {379 errors = (psVector*)psUnaryOp(errors, weights, "sqrt");380 if (params->variances) { 381 errors = (psVector*)psUnaryOp(errors, variances, "sqrt"); 380 382 } 381 383 … … 386 388 outputImage[yOut][xOut] = NAN; 387 389 outputMask[yOut][xOut] = params->blank; 388 if (params-> weights) {389 output Weight[yOut][xOut] = NAN;390 if (params->variances) { 391 outputVariance[yOut][xOut] = NAN; 390 392 } 391 393 sigma->data.F32[yOut][xOut] = NAN; … … 393 395 outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); 394 396 outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank; 395 if (params-> weights) {397 if (params->variances) { 396 398 float stdev = psStatsGetValue(stats, combineStdev); 397 output Weight[yOut][xOut] = PS_SQR(stdev); // Variance399 outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance 398 400 // XXXX this is not the correct formal error. 399 401 // also, the weighted mean is not obviously the correct thing here … … 412 414 psFree(pixels); 413 415 psFree(mask); 414 psFree( weights);416 psFree(variances); 415 417 psFree(errors); 416 418 psFree(stats); -
branches/pap_branch_20090128/psModules/src/imcombine/pmReadoutCombine.h
r21183 r21211 5 5 * @author Paul Price, IfA 6 6 * 7 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $8 * @date $Date: 2009-01-2 7 06:39:38$7 * @version $Revision: 1.15.2.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2009-01-29 00:33:51 $ 9 9 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii 10 10 */ … … 22 22 typedef struct { 23 23 psStatsOptions combine; ///< Statistic to use when performing the combination 24 psImageMaskType maskVal; ///< Mask value25 psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels24 psImageMaskType maskVal; ///< Mask value 25 psImageMaskType blank; ///< Mask value to give blank (i.e., no data) pixels 26 26 int nKeep; ///< Mimimum number of pixels to keep 27 27 float fracHigh; ///< Fraction of high pixels to immediately throw … … 29 29 int iter; ///< Number of iterations for clipping (for CLIPPED_MEAN only) 30 30 float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) 31 bool weights; ///< Use the supplied weights (instead of calculated stdev)?31 bool variances; ///< Use the supplied variances (instead of calculated stdev)? 32 32 } pmCombineParams; 33 33 … … 41 41 /// Combine multiple readouts, applying zero and scale, with optional minmax clipping 42 42 bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned 43 const psArray *inputs, ///< Array of input readouts (F32 image and weight, U8 mask)43 const psArray *inputs, ///< Array of input readouts 44 44 const psVector *zero, ///< Zero corrections to subtract from input, or NULL 45 45 const psVector *scale, ///< Scale corrections to divide into input, or NULL -
branches/pap_branch_20090128/psModules/src/imcombine/pmStack.c
r21183 r21211 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $11 * @date $Date: 2009-01-2 7 06:39:38$10 * @version $Revision: 1.46.2.1 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2009-01-29 00:33:51 $ 12 12 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii 13 13 * … … 263 263 264 264 psImage *image = data->readout->image; // Image of interest 265 psImage *variance = data->readout-> weight; // Variance ("weight")map of interest265 psImage *variance = data->readout->variance; // Variance map of interest 266 266 pixelData->data.F32[num] = image->data.F32[yIn][xIn]; 267 267 if (variance) { … … 442 442 *numCols = data->readout->image->numCols; 443 443 *numRows = data->readout->image->numRows; 444 if (data->readout-> weight) {444 if (data->readout->variance) { 445 445 *haveVariances = true; 446 PS_ASSERT_IMAGE_NON_NULL(data->readout-> weight, false);447 PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout-> weight, false);448 PS_ASSERT_IMAGE_TYPE(data->readout-> weight, PS_TYPE_F32, false);446 PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); 447 PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); 448 PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); 449 449 } 450 450 *haveRejects = (data->reject != NULL); … … 473 473 PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); 474 474 if (*haveVariances) { 475 PS_ASSERT_IMAGE_NON_NULL(data->readout-> weight, false);476 PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout-> weight, false);477 PS_ASSERT_IMAGE_TYPE(data->readout-> weight, PS_TYPE_F32, false);475 PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false); 476 PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false); 477 PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); 478 478 } 479 479 } … … 649 649 psImage *combinedImage = combined->image; // Combined image 650 650 psImage *combinedMask = combined->mask; // Combined mask 651 psImage *combinedVariance = combined-> weight; // Combined variance map651 psImage *combinedVariance = combined->variance; // Combined variance map 652 652 653 653 psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, … … 702 702 } 703 703 704 psImage *combinedVariance = combined-> weight; // Combined variance map704 psImage *combinedVariance = combined->variance; // Combined variance map 705 705 if (haveVariances && !combinedVariance) { 706 combined-> weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);707 combinedVariance = combined-> weight;706 combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 707 combinedVariance = combined->variance; 708 708 } 709 709 -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c
r21183 r21211 60 60 61 61 // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel 62 // This is required to keep the relative scaling between the image and the weightmap62 // This is required to keep the relative scaling between the image and the variance map 63 63 psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); 64 64 … … 287 287 288 288 // Convolve an image using FFT 289 static void convolve WeightFFT(psImage *target,// Place the result in here290 psImage * weight, // Weightmap to convolve289 static void convolveVarianceFFT(psImage *target,// Place the result in here 290 psImage *variance, // Variance map to convolve 291 291 psImage *sys, // Systematic error image 292 292 psImage *mask, // Mask image … … 302 302 bool threaded = pmSubtractionThreaded(); // Are we running threaded? 303 303 304 psImage *sub Weight = convolveSubsetAlloc(weight, border, threaded); // Weightmap304 psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map 305 305 psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image 306 306 psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask 307 307 308 308 // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once 309 psImage *conv Weight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight309 psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance 310 310 psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys 311 311 312 convolveSubsetFree( weight, subWeight, threaded);312 convolveSubsetFree(variance, subVariance, threaded); 313 313 convolveSubsetFree(sys, subSys, threaded); 314 314 convolveSubsetFree(mask, subMask, threaded); … … 319 319 for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { 320 320 for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) { 321 target->data.F32[yTarget][xTarget] = conv Weight->data.F32[ySource][xSource] +321 target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] + 322 322 convSys->data.F32[ySource][xSource]; 323 323 } … … 326 326 int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy 327 327 for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { 328 memcpy(&target->data.F32[yTarget][xMin], &conv Weight->data.F32[ySource][size], numBytes);329 } 330 } 331 332 psFree(conv Weight);328 memcpy(&target->data.F32[yTarget][xMin], &convVariance->data.F32[ySource][size], numBytes); 329 } 330 } 331 332 psFree(convVariance); 333 333 psFree(convSys); 334 334 … … 361 361 // Convolve a region of an image 362 362 static inline void convolveRegion(psImage *convImage, // Convolved image (output) 363 psImage *conv Weight, // Convolved weightmap (output), or NULL363 psImage *convVariance, // Convolved variance map (output), or NULL 364 364 psImage *convMask, // Convolve mask (output), or NULL 365 365 psKernel **kernelImage, // Convolution kernel for the image 366 psKernel **kernel Weight, // Convolution kernel for the weightmap, or NULL366 psKernel **kernelVariance, // Convolution kernel for the variance map, or NULL 367 367 psImage *image, // Image to convolve 368 psImage * weight, // Weightmap to convolve, or NULL368 psImage *variance, // Variance map to convolve, or NULL 369 369 psImage *sys, // Systematic error image, or NULL 370 370 psImage *subMask, // Subtraction mask … … 381 381 { 382 382 *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); 383 if ( weight|| subMask) {384 *kernel Weight = varianceKernel(*kernelWeight, *kernelImage);383 if (variance || subMask) { 384 *kernelVariance = varianceKernel(*kernelVariance, *kernelImage); 385 385 } 386 386 … … 398 398 } 399 399 400 // Convolve the image and weight400 // Convolve the image and variance 401 401 if (useFFT) { 402 402 // Use Fast Fourier Transform to do the convolution 403 403 // This provides a big speed-up for large kernels 404 404 convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); 405 if ( weight) {406 convolve WeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size);405 if (variance) { 406 convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size); 407 407 } 408 408 } else { 409 409 // XXX Direct convolution doesn't account for bad pixels yet 410 410 convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); 411 if ( weight) {412 convolveDirect(conv Weight, weight, *kernelWeight, region, 0.0, kernels->size);411 if (variance) { 412 convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size); 413 413 } 414 414 } … … 885 885 psFree(stamp->image1); 886 886 psFree(stamp->image2); 887 psFree(stamp-> weight);888 stamp->image1 = stamp->image2 = stamp-> weight= NULL;887 psFree(stamp->variance); 888 stamp->image1 = stamp->image2 = stamp->variance = NULL; 889 889 psFree(stamp->matrix1); 890 890 psFree(stamp->matrix2); … … 989 989 990 990 991 // XXX Put kernelImage, kernel Weightand polyValues on thread-dependent data991 // XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data 992 992 static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image 993 993 int x0, int y0, // Offsets for image … … 1010 1010 1011 1011 psKernel *kernelImage = NULL; // Kernel for the images 1012 psKernel *kernel Weight = NULL; // Kernel for the weightmaps1012 psKernel *kernelVariance = NULL; // Kernel for the variance maps 1013 1013 1014 1014 // Only generate polynomial values every kernel footprint, since we have already assumed … … 1020 1020 1021 1021 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1022 convolveRegion(out1->image, out1-> weight, convMask, &kernelImage, &kernelWeight,1023 ro1->image, ro1-> weight, sys1, subMask, kernels, polyValues, background, *region,1022 convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance, 1023 ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region, 1024 1024 maskBad, maskPoor, poorFrac, useFFT, false); 1025 1025 } 1026 1026 if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1027 convolveRegion(out2->image, out2-> weight, convMask, &kernelImage, &kernelWeight,1028 ro2->image, ro2-> weight, sys2, subMask, kernels, polyValues, background, *region,1027 convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance, 1028 ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region, 1029 1029 maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); 1030 1030 } 1031 1031 1032 1032 psFree(kernelImage); 1033 psFree(kernel Weight);1033 psFree(kernelVariance); 1034 1034 psFree(polyValues); 1035 1035 … … 1148 1148 } 1149 1149 } 1150 if (ro1-> weight) {1151 if (!out1-> weight) {1152 out1-> weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);1150 if (ro1->variance) { 1151 if (!out1->variance) { 1152 out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 1153 1153 if (threaded) { 1154 psMutexInit(out1-> weight);1154 psMutexInit(out1->variance); 1155 1155 } 1156 1156 } 1157 psImageInit(out1-> weight, 0.0);1157 psImageInit(out1->variance, 0.0); 1158 1158 } 1159 1159 } … … 1165 1165 } 1166 1166 } 1167 if (ro2-> weight) {1168 if (!out2-> weight) {1169 out2-> weight= psImageAlloc(numCols, numRows, PS_TYPE_F32);1167 if (ro2->variance) { 1168 if (!out2->variance) { 1169 out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 1170 1170 if (threaded) { 1171 psMutexInit(out2-> weight);1171 psMutexInit(out2->variance); 1172 1172 } 1173 1173 } 1174 psImageInit(out2-> weight, 0.0);1174 psImageInit(out2->variance, 0.0); 1175 1175 } 1176 1176 } … … 1232 1232 psImage *polyValues = NULL; // Pre-calculated polynomial values 1233 1233 psKernel *kernelImage = NULL; // Kernel for the images 1234 psKernel *kernel Weight = NULL; // Kernel for the weightmaps1234 psKernel *kernelVariance = NULL; // Kernel for the variance maps 1235 1235 #endif 1236 1236 … … 1345 1345 if (out2) { 1346 1346 out2->image = psMemIncrRefCounter(ro2->image); 1347 out2-> weight = psMemIncrRefCounter(ro2->weight);1347 out2->variance = psMemIncrRefCounter(ro2->variance); 1348 1348 out2->mask = psMemIncrRefCounter(ro2->mask); 1349 1349 } … … 1352 1352 if (out1) { 1353 1353 out1->image = psMemIncrRefCounter(ro1->image); 1354 out1-> weight = psMemIncrRefCounter(ro1->weight);1354 out1->variance = psMemIncrRefCounter(ro1->variance); 1355 1355 out1->mask = psMemIncrRefCounter(ro1->mask); 1356 1356 } -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionEquation.c
r20561 r21211 24 24 static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication 25 25 const psKernel *image2, // Second image in multiplication 26 const psKernel * weight, // Weightimage26 const psKernel *variance, // Variance image 27 27 int footprint // (Half-)Size of stamp 28 28 ) … … 31 31 for (int y = - footprint; y <= footprint; y++) { 32 32 for (int x = - footprint; x <= footprint; x++) { 33 sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x];33 sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; 34 34 } 35 35 } … … 42 42 const psKernel *image1, // First image in multiplication 43 43 const psKernel *image2, // Second image in multiplication 44 const psKernel * weight, // Weightimage44 const psKernel *variance, // Variance image 45 45 const psImage *polyValues, // Spatial polynomial values 46 46 int numKernels, // Number of kernel basis functions … … 50 50 ) 51 51 { 52 double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products52 double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products 53 53 if (!isfinite(sum)) { 54 54 return false; … … 79 79 const psKernel *image1, // First image in multiplication 80 80 const psKernel *image2, // Second image in multiplication 81 const psKernel * weight, // Weightimage81 const psKernel *variance, // Variance image 82 82 const psImage *polyValues, // Spatial polynomial values 83 83 int numKernels, // Number of kernel basis functions … … 87 87 ) 88 88 { 89 double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products89 double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products 90 90 if (!isfinite(sum)) { 91 91 return false; … … 120 120 const psArray *convolutions1, // Convolutions for element 1 121 121 const psArray *convolutions2, // Convolutions for element 2 122 const psKernel * weight, // Weightimage122 const psKernel *variance, // Variance image 123 123 const psImage *polyValues, // Polynomial values 124 124 int numKernels, // Number of kernel basis functions … … 135 135 psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element 136 136 137 if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels,137 if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels, 138 138 footprint, spatialOrder, symmetric)) { 139 139 psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); … … 151 151 const psArray *convolutions, // Convolutions of source with kernels 152 152 const psKernel *input, // Input stamp, or NULL 153 const psKernel * weight, // Weightstamp153 const psKernel *variance, // Variance stamp 154 154 const psImage *polyValues, // Spatial polynomial values 155 155 int footprint, // (Half-)Size of stamp … … 171 171 172 172 // Square part of the matrix (convolution-convolution products) 173 if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels,173 if (!calculateMatrixSquare(matrix, convolutions, convolutions, variance, polyValues, numKernels, 174 174 spatialOrder, footprint)) { 175 175 return false; … … 185 185 186 186 // Normalisation-convolution terms 187 if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels,187 if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, variance, polyValues, numKernels, 188 188 footprint, spatialOrder, true)) { 189 189 psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); … … 195 195 for (int y = - footprint; y <= footprint; y++) { 196 196 for (int x = - footprint; x <= footprint; x++) { 197 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];197 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 198 198 } 199 199 } … … 218 218 for (int y = - footprint; y <= footprint; y++) { 219 219 for (int x = - footprint; x <= footprint; x++) { 220 double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x];220 double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; 221 221 double value = input->kernel[y][x] * invNoise2; 222 222 sumI += value; … … 253 253 const psKernel *input, // Input stamp, or NULL if !normAndBG 254 254 const psKernel *target, // Target stamp 255 const psKernel * weight, // Weightstamp255 const psKernel *variance, // Variance stamp 256 256 const psImage *polyValues, // Spatial polynomial values 257 257 int footprint, // (Half-)Size of stamp … … 277 277 for (int y = - footprint; y <= footprint; y++) { 278 278 for (int x = - footprint; x <= footprint; x++) { 279 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x];279 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 280 280 } 281 281 } … … 297 297 for (int y = - footprint; y <= footprint; y++) { 298 298 for (int x = - footprint; x <= footprint; x++) { 299 float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x];299 float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; 300 300 sumIT += value * input->kernel[y][x]; 301 301 sumT += value; … … 329 329 const psArray *convolutions2, // Convolutions of image 2 330 330 const psKernel *image1, // Image 1 stamp 331 const psKernel * weight, // Weightstamp331 const psKernel *variance, // Variance stamp 332 332 const psImage *polyValues, // Spatial polynomial values 333 333 int footprint // (Half-)Size of stamp … … 348 348 PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); 349 349 350 if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels,350 if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels, 351 351 spatialOrder, footprint)) { 352 352 return false; … … 356 356 // Normalisation 357 357 psKernel *conv = convolutions2->data[i]; // Convolution 358 if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels,358 if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, variance, polyValues, numKernels, 359 359 footprint, spatialOrder, false)) { 360 360 psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); … … 366 366 for (int y = - footprint; y <= footprint; y++) { 367 367 for (int x = - footprint; x <= footprint; x++) { 368 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];368 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 369 369 } 370 370 } … … 559 559 case PM_SUBTRACTION_MODE_1: 560 560 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 561 stamp-> weight, polyValues, footprint, true);561 stamp->variance, polyValues, footprint, true); 562 562 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 563 stamp->image2, stamp-> weight, polyValues, footprint, true);563 stamp->image2, stamp->variance, polyValues, footprint, true); 564 564 break; 565 565 case PM_SUBTRACTION_MODE_2: 566 566 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, 567 stamp-> weight, polyValues, footprint, true);567 stamp->variance, polyValues, footprint, true); 568 568 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, 569 stamp->image1, stamp-> weight, polyValues, footprint, true);569 stamp->image1, stamp->variance, polyValues, footprint, true); 570 570 break; 571 571 case PM_SUBTRACTION_MODE_DUAL: … … 581 581 #endif 582 582 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 583 stamp-> weight, polyValues, footprint, true);583 stamp->variance, polyValues, footprint, true); 584 584 status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, 585 stamp-> weight, polyValues, footprint, false);585 stamp->variance, polyValues, footprint, false); 586 586 status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, 587 stamp->convolutions2, stamp->image1, stamp-> weight, polyValues,587 stamp->convolutions2, stamp->image1, stamp->variance, polyValues, 588 588 footprint); 589 589 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 590 stamp->image2, stamp-> weight, polyValues, footprint, true);590 stamp->image2, stamp->variance, polyValues, footprint, true); 591 591 status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, 592 stamp->image2, stamp-> weight, polyValues, footprint, false);592 stamp->image2, stamp->variance, polyValues, footprint, false); 593 593 break; 594 594 default: … … 1033 1033 1034 1034 // Calculate residuals 1035 psKernel * weight = stamp->weight; // Weightpostage stamp1035 psKernel *variance = stamp->variance; // Variance postage stamp 1036 1036 psImageInit(residual->image, 0.0); 1037 1037 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { … … 1098 1098 for (int y = - footprint; y <= footprint; y++) { 1099 1099 for (int x = - footprint; x <= footprint; x++) { 1100 double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x];1100 double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x]; 1101 1101 deviation += dev; 1102 1102 #ifdef TESTING … … 1141 1141 psFitsClose(fits); 1142 1142 } 1143 if (stamp-> weight) {1143 if (stamp->variance) { 1144 1144 psString filename = NULL; 1145 psStringAppend(&filename, "stamp_ weight_%03d.fits", i);1145 psStringAppend(&filename, "stamp_variance_%03d.fits", i); 1146 1146 psFits *fits = psFitsOpen(filename, "w"); 1147 1147 psFree(filename); 1148 psFitsWriteImage(fits, NULL, stamp-> weight->image, 0, NULL);1148 psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL); 1149 1149 psFitsClose(fits); 1150 1150 } -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionMatch.c
r21183 r21211 60 60 const pmReadout *ro2, // Readout 2 61 61 const psImage *subMask, // Mask for subtraction, or NULL 62 psImage * weight, // Weightmap62 psImage *variance, // Variance map 63 63 const psRegion *region, // Region of interest, or NULL 64 64 float threshold, // Threshold for stamp finding … … 80 80 81 81 psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); 82 if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) {82 if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) { 83 83 psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); 84 84 return false; … … 111 111 conv1->mask = NULL; 112 112 } 113 if (conv1-> weight) {114 psFree(conv1-> weight);115 conv1-> weight= NULL;113 if (conv1->variance) { 114 psFree(conv1->variance); 115 conv1->variance = NULL; 116 116 } 117 117 } … … 126 126 conv2->mask = NULL; 127 127 } 128 if (conv2-> weight) {129 psFree(conv2-> weight);130 conv2-> weight= NULL;128 if (conv2->variance) { 129 psFree(conv2->variance); 130 conv2->variance = NULL; 131 131 } 132 132 } … … 191 191 } 192 192 193 // Where does our weightmap come from?194 // Getting the weightexactly right is not necessary --- it's just used for weighting.195 psImage * weight = NULL; // Weightimage to use196 if (ro1-> weight && ro2->weight) {197 weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight);198 } else if (ro1-> weight) {199 weight = psMemIncrRefCounter(ro1->weight);200 } else if (ro2-> weight) {201 weight = psMemIncrRefCounter(ro2->weight);193 // Where does our variance map come from? 194 // Getting the variance exactly right is not necessary --- it's just used for weighting. 195 psImage *variance = NULL; // Variance image to use 196 if (ro1->variance && ro2->variance) { 197 variance = (psImage*)psBinaryOp(NULL, ro1->variance, "+", ro2->variance); 198 } else if (ro1->variance) { 199 variance = psMemIncrRefCounter(ro1->variance); 200 } else if (ro2->variance) { 201 variance = psMemIncrRefCounter(ro2->variance); 202 202 } else { 203 weight= (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);203 variance = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); 204 204 } 205 205 … … 274 274 // We get the stamps here; we will also attempt to get stamps at the first iteration, but it 275 275 // doesn't matter. 276 if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,276 if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, threshold, stampSpacing, 277 277 size, footprint, subMode)) { 278 278 goto MATCH_ERROR; … … 336 336 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); 337 337 338 if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,338 if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing, 339 339 size, footprint, subMode)) { 340 340 goto MATCH_ERROR; … … 439 439 psFree(subMask); 440 440 subMask = NULL; 441 psFree( weight);442 weight= NULL;443 444 if (!pmSubtractionBorder(conv1->image, conv1-> weight, conv1->mask, size, maskBad)) {441 psFree(variance); 442 variance = NULL; 443 444 if (!pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) { 445 445 psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); 446 446 goto MATCH_ERROR; … … 482 482 psFree(kernels); 483 483 psFree(stamps); 484 psFree( weight);484 psFree(variance); 485 485 psFree(rng); 486 486 return false; -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionParams.c
r18287 r21211 71 71 double *sumII, // Sum of I(x)^2/sigma(x)^2 72 72 double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 73 const pmSubtractionStamp *stamp, // Stamp with weight73 const pmSubtractionStamp *stamp, // Stamp with variance 74 74 const psKernel *target, // Target stamp 75 75 int kernelIndex, // Index for kernel component … … 78 78 ) 79 79 { 80 psKernel * weight = stamp->weight; // Weight, sigma(x)^280 psKernel *variance = stamp->variance; // Variance, sigma(x)^2 81 81 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 82 82 83 83 for (int y = -footprint; y <= footprint; y++) { 84 84 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 85 psF32 *wt = & weight->kernel[y][-footprint]; // Dereference weight85 psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance 86 86 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 87 87 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { … … 98 98 static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 99 99 double *sumCC, // Sum of conv(x)^2/sigma(x)^2 100 const pmSubtractionStamp *stamp, // Stamp with input and weight100 const pmSubtractionStamp *stamp, // Stamp with input and variance 101 101 int kernelIndex, // Index for kernel component 102 102 int footprint, // Size of region of interest … … 104 104 ) 105 105 { 106 psKernel * weight = stamp->weight; // Weight, sigma(x)^2106 psKernel *variance = stamp->variance; // Variance, sigma(x)^2 107 107 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 108 108 109 109 for (int y = -footprint; y <= footprint; y++) { 110 psF32 *wt = & weight->kernel[y][-footprint]; // Dereference weight110 psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance 111 111 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 112 112 for (int x = -footprint; x <= footprint; x++, wt++, conv++) { … … 120 120 121 121 static double accumulateChi2(const psKernel *target, // Target stamp 122 pmSubtractionStamp *stamp, // Stamp with weight122 pmSubtractionStamp *stamp, // Stamp with variance 123 123 int kernelIndex, // Index for kernel component 124 124 double coeff, // Coefficient of convolution … … 129 129 { 130 130 double chi2 = 0.0; 131 psKernel * weight = stamp->weight; // Weight, sigma(x)^2131 psKernel *variance = stamp->variance; // Variance, sigma(x)^2 132 132 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 133 133 134 134 for (int y = -footprint; y <= footprint; y++) { 135 135 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 136 psF32 *wt = & weight->kernel[y][-footprint]; // Dereference weight136 psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance 137 137 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 138 138 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { … … 146 146 // Return the initial value of chi^2 147 147 static double initialChi2(const psKernel *target, // Target stamp 148 const pmSubtractionStamp *stamp, // Stamp with weight148 const pmSubtractionStamp *stamp, // Stamp with variance 149 149 int footprint, // Size of convolution 150 150 pmSubtractionMode mode // Mode of subtraction 151 151 ) 152 152 { 153 psKernel * weight = stamp->weight; // Weightmap153 psKernel *variance = stamp->variance; // Variance map 154 154 psKernel *source; // Source stamp 155 155 switch (mode) { … … 167 167 for (int y = -footprint; y <= footprint; y++) { 168 168 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 169 psF32 *wt = & weight->kernel[y][-footprint]; // Dereference weight169 psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance 170 170 psF32 *ref = &source->kernel[y][-footprint]; // Derference reference 171 171 for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { … … 180 180 // Subtract a convolution from the input 181 181 static void subtractConvolution(psKernel *target, // Target stamp 182 const pmSubtractionStamp *stamp, // Stamp with weight182 const pmSubtractionStamp *stamp, // Stamp with variance 183 183 int kernelIndex, // Index for kernel component 184 184 float coeff, // Coefficient of subtraction … … 288 288 289 289 // This sum is invariant to the kernel 290 psKernel * weight = stamp->weight; // Weightmap for stamp290 psKernel *variance = stamp->variance; // Variance map for stamp 291 291 for (int v = -footprint; v <= footprint; v++) { 292 psF32 *wt = & weight->kernel[v][-footprint]; // Dereference weightmap292 psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map 293 293 for (int u = -footprint; u <= footprint; u++, wt++) { 294 294 sum1 += 1.0 / *wt; -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionStamps.c
r21183 r21211 54 54 psFree(stamp->image1); 55 55 psFree(stamp->image2); 56 psFree(stamp-> weight);56 psFree(stamp->variance); 57 57 psFree(stamp->convolutions1); 58 58 psFree(stamp->convolutions2); … … 211 211 stamp->image1 = NULL; 212 212 stamp->image2 = NULL; 213 stamp-> weight= NULL;213 stamp->variance = NULL; 214 214 stamp->convolutions1 = NULL; 215 215 stamp->convolutions2 = NULL; … … 333 333 psFree(stamp->image1); 334 334 psFree(stamp->image2); 335 psFree(stamp-> weight);335 psFree(stamp->variance); 336 336 psFree(stamp->convolutions1); 337 337 psFree(stamp->convolutions2); 338 stamp->image1 = stamp->image2 = stamp-> weight= NULL;338 stamp->image1 = stamp->image2 = stamp->variance = NULL; 339 339 stamp->convolutions1 = stamp->convolutions2 = NULL; 340 340 … … 480 480 481 481 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, 482 psImage * weight, int kernelSize)482 psImage *variance, int kernelSize) 483 483 { 484 484 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); … … 490 490 PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); 491 491 } 492 PS_ASSERT_IMAGE_NON_NULL( weight, false);493 PS_ASSERT_IMAGES_SIZE_EQUAL( weight, image1, false);494 PS_ASSERT_IMAGE_TYPE( weight, PS_TYPE_F32, false);492 PS_ASSERT_IMAGE_NON_NULL(variance, false); 493 PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); 494 PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); 495 495 PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); 496 496 … … 520 520 assert(stamp->image1 == NULL); 521 521 assert(stamp->image2 == NULL); 522 assert(stamp-> weight== NULL);522 assert(stamp->variance == NULL); 523 523 524 524 psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest … … 534 534 } 535 535 536 psImage *wtSub = psImageSubset( weight, region); // Subimage with stamp537 stamp-> weight= psKernelAllocFromImage(wtSub, size, size);536 psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp 537 stamp->variance = psKernelAllocFromImage(wtSub, size, size); 538 538 psFree(wtSub); // Drop reference 539 539 -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionStamps.h
r20465 r21211 60 60 psKernel *image1; ///< Reference image postage stamp 61 61 psKernel *image2; ///< Input image postage stamp 62 psKernel * weight; ///< Weightimage postage stamp, or NULL62 psKernel *variance; ///< Variance image postage stamp, or NULL 63 63 psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL 64 64 psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL … … 120 120 psImage *image1, ///< Reference image 121 121 psImage *image2, ///< Input image (or NULL) 122 psImage * weight, ///< Weight (variance)map122 psImage *variance, ///< Variance map 123 123 int kernelSize ///< Kernel half-size 124 124 ); -
branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionThreads.c
r19765 r21211 27 27 psMutexInit(ro->image); 28 28 } 29 if (ro-> weight) {30 psMutexInit(ro-> weight);29 if (ro->variance) { 30 psMutexInit(ro->variance); 31 31 } 32 32 … … 43 43 psMutexDestroy(ro->image); 44 44 } 45 if (ro-> weight) {46 psMutexDestroy(ro-> weight);45 if (ro->variance) { 46 psMutexDestroy(ro->variance); 47 47 } 48 48 -
branches/pap_branch_20090128/psModules/src/objects/pmSource.c
r21183 r21211 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.67 $ $Name: not supported by cvs2svn $9 * @date $Date: 2009-01-2 7 06:39:38$8 * @version $Revision: 1.67.2.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 44 44 psFree(tmp->peak); 45 45 psFree(tmp->pixels); 46 psFree(tmp-> weight);46 psFree(tmp->variance); 47 47 psFree(tmp->maskObj); 48 48 psFree(tmp->maskView); … … 66 66 67 67 psFree (source->pixels); 68 psFree (source-> weight);68 psFree (source->variance); 69 69 psFree (source->maskObj); 70 70 psFree (source->maskView); … … 73 73 74 74 source->pixels = NULL; 75 source-> weight= NULL;75 source->variance = NULL; 76 76 source->maskObj = NULL; 77 77 source->maskView = NULL; … … 101 101 source->peak = NULL; 102 102 source->pixels = NULL; 103 source-> weight= NULL;103 source->variance = NULL; 104 104 source->maskObj = NULL; 105 105 source->maskView = NULL; … … 144 144 } 145 145 // this copy is used to allow multiple fit attempts on the 146 // same object. the pixels, weight, and mask arrays all point to146 // same object. the pixels, variance, and mask arrays all point to 147 147 // the same original subarrays. the peak and moments point at 148 148 // the original values. … … 167 167 // We want a new view, but pointing at the same pixels. 168 168 source->pixels = psImageCopyView(NULL, in->pixels); 169 source-> weight = psImageCopyView(NULL, in->weight);169 source->variance = psImageCopyView(NULL, in->variance); 170 170 source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL; 171 171 … … 199 199 // these images are subset images of the equivalent parents 200 200 mySource->pixels = psImageSubset(readout->image, srcRegion); 201 if (readout-> weight) {202 mySource-> weight = psImageSubset(readout->weight, srcRegion);201 if (readout->variance) { 202 mySource->variance = psImageSubset(readout->variance, srcRegion); 203 203 } 204 204 if (readout->mask) { … … 237 237 238 238 extend |= (mySource->pixels == NULL); 239 extend |= (mySource-> weight== NULL);239 extend |= (mySource->variance == NULL); 240 240 extend |= (mySource->maskObj == NULL); 241 241 extend |= (mySource->maskView == NULL); … … 245 245 // re-create the subimage 246 246 psFree (mySource->pixels); 247 psFree (mySource-> weight);247 psFree (mySource->variance); 248 248 psFree (mySource->maskView); 249 249 250 250 mySource->pixels = psImageSubset(readout->image, newRegion); 251 mySource-> weight = psImageSubset(readout->weight, newRegion);251 mySource->variance = psImageSubset(readout->variance, newRegion); 252 252 mySource->maskView = psImageSubset(readout->mask, newRegion); 253 253 mySource->region = newRegion; … … 670 670 pmSource->peak 671 671 pmSource->pixels 672 pmSource-> weight672 pmSource->variance 673 673 pmSource->mask 674 674 … … 737 737 738 738 psF32 *vPix = source->pixels->data.F32[row]; 739 psF32 *vWgt = source-> weight->data.F32[row];739 psF32 *vWgt = source->variance->data.F32[row]; 740 740 psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; 741 741 … … 923 923 if (mode & PM_MODEL_OP_NOISE) { 924 924 // XXX need to scale by the gain... 925 target = source-> weight->data.F32;925 target = source->variance->data.F32; 926 926 } 927 927 … … 945 945 psImage *target = source->pixels; 946 946 if (mode & PM_MODEL_OP_NOISE) { 947 target = source-> weight;947 target = source->variance; 948 948 } 949 949 -
branches/pap_branch_20090128/psModules/src/objects/pmSource.h
r21183 r21211 3 3 * @author EAM, IfA; GLG, MHPCC 4 4 * 5 * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $6 * @date $Date: 2009-01-2 7 06:39:38$5 * @version $Revision: 1.27.2.1 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2009-01-29 00:33:51 $ 7 7 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 8 8 */ … … 72 72 pmPeak *peak; ///< Description of peak pixel. 73 73 psImage *pixels; ///< Rectangular region including object pixels. 74 psImage * weight; ///< Image variance.74 psImage *variance; ///< Image variance. 75 75 psImage *maskObj; ///< unique mask for this object which marks included pixels associated with objects. 76 76 psImage *maskView; ///< view into global image mask for this object region … … 210 210 psMetadata *metadata, ///< Contains classification parameters 211 211 pmPSFClump clump, ///< Statistics about the PSF clump 212 psImageMaskType maskSat ///< Mask value for saturated pixels212 psImageMaskType maskSat ///< Mask value for saturated pixels 213 213 ); 214 214 -
branches/pap_branch_20090128/psModules/src/objects/pmSourceFitModel.c
r21183 r21211 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $9 * @date $Date: 2009-01-2 7 06:39:38$8 * @version $Revision: 1.29.2.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 62 62 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 63 63 PS_ASSERT_PTR_NON_NULL(source->maskObj, false); 64 PS_ASSERT_PTR_NON_NULL(source-> weight, false);64 PS_ASSERT_PTR_NON_NULL(source->variance, false); 65 65 66 66 psBool fitStatus = true; … … 84 84 continue; 85 85 } 86 // skip zero- weightpoints87 if (source-> weight->data.F32[i][j] == 0) {86 // skip zero-variance points 87 if (source->variance->data.F32[i][j] == 0) { 88 88 continue; 89 89 } … … 102 102 103 103 // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky 104 // as weightto avoid the bias from systematic errors here we would just use the104 // as variance to avoid the bias from systematic errors here we would just use the 105 105 // source sky variance 106 106 if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { 107 yErr->data.F32[nPix] = 1.0 / source-> weight->data.F32[i][j];107 yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; 108 108 } else { 109 109 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; -
branches/pap_branch_20090128/psModules/src/objects/pmSourceFitSet.c
r21183 r21211 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $9 * @date $Date: 2009-01-2 7 06:39:38$8 * @version $Revision: 1.14.2.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 11 11 * … … 55 55 56 56 if (!fitSets) { 57 fitSets = psArrayAlloc (PS_MAX (1, nThreads));58 } 59 60 // the allocated elements should be NULL on psArrayAlloc, 61 // and a previously allocated array of fitSets should have been cleared 57 fitSets = psArrayAlloc (PS_MAX (1, nThreads)); 58 } 59 60 // the allocated elements should be NULL on psArrayAlloc, 61 // and a previously allocated array of fitSets should have been cleared 62 62 // before pmSourceFitSetInit is called 63 63 for (int i = 0; i < fitSets->n; i++) { 64 psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?");64 psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?"); 65 65 } 66 66 return true; … … 125 125 // we do not need to lock to do this.... 126 126 for (int i = 0; i < fitSets->n; i++) { 127 thisSet = fitSets->data[i];128 if (!thisSet) continue;129 if (thisSet->thread == id) {130 break;131 }132 thisSet = NULL;127 thisSet = fitSets->data[i]; 128 if (!thisSet) continue; 129 if (thisSet->thread == id) { 130 break; 131 } 132 thisSet = NULL; 133 133 } 134 134 psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared"); … … 137 137 138 138 pthread_mutex_lock (&fitSetInitMutex); 139 139 140 140 // find an entry that is NULL and place it there 141 141 for (int i = 0; i < fitSets->n; i++) { 142 if (fitSets->data[i]) continue;143 fitSets->data[i] = thisSet;144 pthread_mutex_unlock (&fitSetInitMutex);145 return thisSet;142 if (fitSets->data[i]) continue; 143 fitSets->data[i] = thisSet; 144 pthread_mutex_unlock (&fitSetInitMutex); 145 return thisSet; 146 146 } 147 147 pthread_mutex_unlock (&fitSetInitMutex); … … 160 160 pmSourceFitSetData *thisSet = NULL; 161 161 for (int i = 0; i < fitSets->n; i++) { 162 thisSet = fitSets->data[i];163 if (!thisSet) continue;164 if (thisSet->thread == id) {165 break;166 }167 thisSet = NULL;162 thisSet = fitSets->data[i]; 163 if (!thisSet) continue; 164 if (thisSet->thread == id) { 165 break; 166 } 167 thisSet = NULL; 168 168 } 169 169 psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found"); … … 184 184 pmSourceFitSetData *thisSet = NULL; 185 185 for (i = 0; i < fitSets->n; i++) { 186 thisSet = fitSets->data[i];187 if (!thisSet) continue;188 if (thisSet->thread == id) {189 break;190 }191 thisSet = NULL;186 thisSet = fitSets->data[i]; 187 if (!thisSet) continue; 188 if (thisSet->thread == id) { 189 break; 190 } 191 thisSet = NULL; 192 192 } 193 193 psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found"); … … 256 256 psVector *derivOne = set->derivSet->data[i]; 257 257 258 // one or the other (param or deriv) must be set259 assert ((deriv != NULL) || (param != NULL));260 261 // if we are setting derive, derivOne and paramOne must be same length258 // one or the other (param or deriv) must be set 259 assert ((deriv != NULL) || (param != NULL)); 260 261 // if we are setting derive, derivOne and paramOne must be same length 262 262 assert ((deriv == NULL) || (paramOne->n == derivOne->n)); 263 263 264 264 for (int j = 0; j < paramOne->n; j++, n++) { 265 if (param) {266 param->data.F32[n] = paramOne->data.F32[j];267 }265 if (param) { 266 param->data.F32[n] = paramOne->data.F32[j]; 267 } 268 268 if (deriv) { 269 269 deriv->data.F32[n] = derivOne->data.F32[j]; … … 275 275 276 276 // distribute parameters from single param and deriv vectors into FitSet models 277 bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) 277 bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) 278 278 { 279 279 PS_ASSERT_VECTOR_NON_NULL(param, false); … … 301 301 302 302 // set the model parameters for this fit set 303 bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, 303 bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, 304 304 const psVector *param, pmSource *source, 305 305 psMinimization *myMin, int nPix, bool fitStatus) … … 420 420 421 421 for (int i = 0; i < modelSet->n; i++) { 422 pmModel *model = modelSet->data[i];423 int nParams = pmModelClassParameterCount (model->type);424 for (int j = 0; j < nParams; j++) {425 fprintf (file, "%d %d : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]);426 }422 pmModel *model = modelSet->data[i]; 423 int nParams = pmModelClassParameterCount (model->type); 424 for (int j = 0; j < nParams; j++) { 425 fprintf (file, "%d %d : %f %f\n", i, j, model->params->data.F32[j], model->dparams->data.F32[j]); 426 } 427 427 } 428 428 return true; … … 435 435 psVector *derivOne = set->derivSet->data[i]; 436 436 for (int j = 0; j < paramOne->n; j++) { 437 fprintf (file, "%d %d : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]);437 fprintf (file, "%d %d : %f %f\n", i, j, paramOne->data.F32[j], derivOne->data.F32[j]); 438 438 } 439 439 } … … 450 450 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 451 451 PS_ASSERT_PTR_NON_NULL(source->maskObj, false); 452 PS_ASSERT_PTR_NON_NULL(source-> weight, false);452 PS_ASSERT_PTR_NON_NULL(source->variance, false); 453 453 454 454 bool fitStatus = true; … … 471 471 continue; 472 472 } 473 // skip zero- weightpoints474 if (source-> weight->data.F32[i][j] == 0) {473 // skip zero-variance points 474 if (source->variance->data.F32[i][j] == 0) { 475 475 continue; 476 476 } … … 489 489 490 490 // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky 491 // as weightto avoid the bias from systematic errors here we would just use the491 // as variance to avoid the bias from systematic errors here we would just use the 492 492 // source sky variance 493 493 if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { 494 yErr->data.F32[nPix] = 1.0 / source-> weight->data.F32[i][j];494 yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; 495 495 } else { 496 496 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; … … 543 543 psFree (params); 544 544 psFree(constraint); 545 pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets545 pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets 546 546 return(false); 547 547 } -
branches/pap_branch_20090128/psModules/src/objects/pmSourceMoments.c
r21183 r21211 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $9 * @date $Date: 2009-01-2 7 06:39:38$8 * @version $Revision: 1.7.2.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 47 47 pmSource->peak 48 48 pmSource->pixels 49 pmSource-> weight(optional)49 pmSource->variance (optional) 50 50 pmSource->mask (optional) 51 51 … … 99 99 100 100 psF32 *vPix = source->pixels->data.F32[row]; 101 psF32 *vWgt = source-> weight->data.F32[row];101 psF32 *vWgt = source->variance->data.F32[row]; 102 102 psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; 103 103 … … 110 110 vMsk++; 111 111 } 112 if (isnan(*vPix)) continue;112 if (isnan(*vPix)) continue; 113 113 114 114 psF32 xDiff = col - xPeak; … … 189 189 190 190 psF32 *vPix = source->pixels->data.F32[row]; 191 psF32 *vWgt = source-> weight->data.F32[row];191 psF32 *vWgt = source->variance->data.F32[row]; 192 192 psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; 193 193 … … 200 200 vMsk++; 201 201 } 202 if (isnan(*vPix)) continue;202 if (isnan(*vPix)) continue; 203 203 204 204 psF32 xDiff = col - xCM; … … 206 206 207 207 // radius is just a function of (xDiff, yDiff) 208 psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);209 psF32 r = sqrt(r2);208 psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); 209 psF32 r = sqrt(r2); 210 210 if (r > radius) continue; 211 211 … … 278 278 279 279 psTrace ("psModules.objects", 4, "Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", 280 source->moments->Mxx, source->moments->Mxy, source->moments->Myy, 280 source->moments->Mxx, source->moments->Mxy, source->moments->Myy, 281 281 source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, 282 282 source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); … … 286 286 287 287 if (source->moments->Mxx < 0) { 288 fprintf (stderr, "error: neg second moment??\n");288 fprintf (stderr, "error: neg second moment??\n"); 289 289 } 290 290 if (source->moments->Myy < 0) { 291 fprintf (stderr, "error: neg second moment??\n");291 fprintf (stderr, "error: neg second moment??\n"); 292 292 } 293 293 … … 341 341 342 342 psF32 *vPix = source->pixels->data.F32[row]; 343 psF32 *vWgt = source-> weight->data.F32[row];343 psF32 *vWgt = source->variance->data.F32[row]; 344 344 psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; 345 345 … … 352 352 vMsk++; 353 353 } 354 if (isnan(*vPix)) continue;354 if (isnan(*vPix)) continue; 355 355 356 356 psF32 xDiff = col - xPeak; … … 407 407 } 408 408 409 # if (PS_TRACE_ON) 409 # if (PS_TRACE_ON) 410 410 float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0); 411 411 float Sxy = XY / Sum; -
branches/pap_branch_20090128/psModules/src/objects/pmSourcePhotometry.c
r21183 r21211 3 3 * @author EAM, IfA; GLG, MHPCC 4 4 * 5 * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $6 * @date $Date: 2009-01-2 7 06:39:38$5 * @version $Revision: 1.48.2.1 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2009-01-29 00:33:51 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 446 446 assert (Pj != NULL); 447 447 448 const psImage *Wi = Mi-> weight;448 const psImage *Wi = Mi->variance; 449 449 if (!unweighted_sum) { 450 450 assert (Wi != NULL); … … 511 511 assert (Pj != NULL); 512 512 513 const psImage *Wi = Mi-> weight;513 const psImage *Wi = Mi->variance; 514 514 if (!unweighted_sum) { 515 515 assert (Wi != NULL); … … 567 567 const psImage *Pi = Mi->pixels; 568 568 assert (Pi != NULL); 569 const psImage *Wi = Mi-> weight;569 const psImage *Wi = Mi->variance; 570 570 if (!unweighted_sum) { 571 571 assert (Wi != NULL); … … 612 612 # endif 613 613 614 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage * weight,614 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, 615 615 psImageMaskType maskVal) 616 616 { … … 618 618 PS_ASSERT_PTR_NON_NULL(image, false); 619 619 PS_ASSERT_PTR_NON_NULL(mask, false); 620 PS_ASSERT_PTR_NON_NULL( weight, false);620 PS_ASSERT_PTR_NON_NULL(variance, false); 621 621 622 622 double dC = 0.0; … … 626 626 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal) 627 627 continue; 628 if ( weight->data.F32[j][i] <= 0)629 continue; 630 dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];628 if (variance->data.F32[j][i] <= 0) 629 continue; 630 dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i]; 631 631 Npix ++; 632 632 } … … 648 648 const psImage *Pi = Mi->modelFlux; 649 649 assert (Pi != NULL); 650 const psImage *Wi = Mi-> weight;650 const psImage *Wi = Mi->variance; 651 651 if (!unweighted_sum) { 652 652 assert (Wi != NULL); … … 706 706 assert (Pj != NULL); 707 707 708 const psImage *Wi = Mi-> weight;708 const psImage *Wi = Mi->variance; 709 709 if (!unweighted_sum) { 710 710 assert (Wi != NULL); … … 770 770 assert (Pj != NULL); 771 771 772 const psImage *Wi = Mi-> weight;772 const psImage *Wi = Mi->variance; 773 773 if (!unweighted_sum) { 774 774 assert (Wi != NULL); -
branches/pap_branch_20090128/psModules/src/objects/pmSourceSky.c
r21183 r21211 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $9 * @date $Date: 2009-01-2 7 06:39:38$8 * @version $Revision: 1.19.2.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2009-01-29 00:33:51 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 113 113 psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__); 114 114 PS_ASSERT_PTR_NON_NULL(source, false); 115 PS_ASSERT_IMAGE_NON_NULL(source-> weight, false);115 PS_ASSERT_IMAGE_NON_NULL(source->variance, false); 116 116 PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false); 117 117 PS_ASSERT_PTR_NON_NULL(source->peak, false); … … 127 127 } 128 128 129 psImage *image = source-> weight;129 psImage *image = source->variance; 130 130 psImage *mask = source->maskObj; 131 131 pmPeak *peak = source->peak;
Note:
See TracChangeset
for help on using the changeset viewer.
