IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21211


Ignore:
Timestamp:
Jan 28, 2009, 2:33:51 PM (17 years ago)
Author:
Paul Price
Message:

Changing pmReadout.weight to variance. Adding pmReadout.covariance which will carry around a covariance pseudo-matrix, allowing us to calculate the pixel-to-pixel variance. pmReadout.covariance now exists and is set, but no mechanism yet to read/write, or use it.

Location:
branches/pap_branch_20090128/psModules/src
Files:
41 edited

Legend:

Unmodified
Added
Removed
  • branches/pap_branch_20090128/psModules/src/camera/pmFPA.c

    r19013 r21211  
    3232    psFree(readout->image);
    3333    psFree(readout->mask);
    34     psFree(readout->weight);
     34    psFree(readout->variance);
     35    psFree(readout->covariance);
    3536    psFree(readout->analysis);
    3637    psFree(readout->bias);
     
    163164    psFree(readout->image);
    164165    psFree(readout->mask);
    165     psFree(readout->weight);
     166    psFree(readout->variance);
     167    psFree(readout->covariance);
    166168    psFree(readout->bias);
    167169
     
    169171
    170172    readout->image = NULL;
    171     readout->weight = NULL;
     173    readout->variance = NULL;
     174    readout->covariance = NULL;
    172175    readout->mask = NULL;
    173176
     
    179182    readout->thisImageScan = 0;
    180183    readout->thisMaskScan = 0;
    181     readout->thisWeightScan = 0;
     184    readout->thisVarianceScan = 0;
    182185
    183186    readout->lastImageScan = 0;
    184187    readout->lastMaskScan = 0;
    185     readout->lastWeightScan = 0;
     188    readout->lastVarianceScan = 0;
    186189}
    187190
     
    197200    if (cell->hdu) {
    198201        psFree(cell->hdu->images);
    199         psFree(cell->hdu->weights);
     202        psFree(cell->hdu->variances);
    200203        psFree(cell->hdu->masks);
    201204        // psFree(cell->hdu->header);
    202205
    203206        cell->hdu->images = NULL;
    204         cell->hdu->weights = NULL;
     207        cell->hdu->variances = NULL;
    205208        cell->hdu->masks = NULL;
    206209        // cell->hdu->header = NULL;
     
    219222    if (chip->hdu) {
    220223        psFree(chip->hdu->images);
    221         psFree(chip->hdu->weights);
     224        psFree(chip->hdu->variances);
    222225        psFree(chip->hdu->masks);
    223226        // psFree(chip->hdu->header);
    224227
    225228        chip->hdu->images = NULL;
    226         chip->hdu->weights = NULL;
     229        chip->hdu->variances = NULL;
    227230        chip->hdu->masks = NULL;
    228231        // chip->hdu->header = NULL;
     
    241244    if (fpa->hdu) {
    242245        psFree(fpa->hdu->images);
    243         psFree(fpa->hdu->weights);
     246        psFree(fpa->hdu->variances);
    244247        psFree(fpa->hdu->masks);
    245248        // psFree(fpa->hdu->header);
    246249
    247250        fpa->hdu->images = NULL;
    248         fpa->hdu->weights = NULL;
     251        fpa->hdu->variances = NULL;
    249252        fpa->hdu->masks = NULL;
    250253        // fpa->hdu->header = NULL;
     
    259262    tmpReadout->image = NULL;
    260263    tmpReadout->mask = NULL;
    261     tmpReadout->weight = NULL;
     264    tmpReadout->variance = NULL;
     265    tmpReadout->covariance = NULL;
    262266    tmpReadout->bias = psListAlloc(NULL);
    263267    tmpReadout->analysis = psMetadataAlloc();
     
    276280    tmpReadout->thisImageScan = 0;
    277281    tmpReadout->thisMaskScan = 0;
    278     tmpReadout->thisWeightScan = 0;
     282    tmpReadout->thisVarianceScan = 0;
    279283
    280284    tmpReadout->lastImageScan = 0;
    281285    tmpReadout->lastMaskScan = 0;
    282     tmpReadout->lastWeightScan = 0;
     286    tmpReadout->lastVarianceScan = 0;
    283287
    284288    tmpReadout->forceScan = false;
  • branches/pap_branch_20090128/psModules/src/camera/pmFPA.h

    r19013 r21211  
    66 * @author Eugene Magnier, IfA
    77 *
    8  * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2008-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 $
    1010 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    1111 */
     
    113113///
    114114/// 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 weight
     115/// or one of multiple coadds).  It contains the actual pixels used in analysis (along with mask and variance
    116116/// maps).  When reading from a FITS file, the images are subimages (from CELL.TRIMSEC) of the pixels read
    117117/// from the appropriate HDU (at the FPA, chip or cell level).  The readout also contains a list of bias
     
    123123    psImage *image;                     ///< Imaging area of readout (corresponds to CELL.TRIMSEC region)
    124124    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)
    126127    psList *bias;                       ///< List of bias (prescan/overscan) images
    127128    psMetadata *analysis;               ///< Readout-level analysis metadata
     
    130131    bool file_exists;                   ///< Does the file for this readout exist (read case only)?
    131132    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 thisWeightScan;                 ///< start scan for next/current read
    137     int lastWeightScan;                 ///< start scan of the last read
     133    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
    138139    bool forceScan;                     ///< Force pmFPARead to obey the above commanded this and last scans.
    139140} pmReadout;
     
    217218        } \
    218219    } \
    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); \
    225228            return RETVAL; \
    226229        } \
     
    236239    PS_ASSERT_IMAGE_NON_NULL((READOUT)->mask, RETVAL);
    237240
    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);
    241248
    242249/// @}
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAConstruct.c

    r17911 r21211  
    14991499                psImage *image = readout->image; // The image
    15001500                psImage *mask = readout->mask; // The mask
    1501                 psImage *weight = readout->weight; // The weight
     1501                psImage *variance = readout->variance; // The variance
    15021502                psList *bias = readout->bias; // The list of bias images
    15031503                if (image) {
    15041504                    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);
    15081509                }
    15091510                if (bias) {
     
    15121513                    while ((biasImage = psListGetAndIncrement(biasIter))) {
    15131514                        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);
    15171519                    }
    15181520                    psFree(biasIter);
     
    15201522                if (mask) {
    15211523                    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);
    15251528                }
    1526                 if (weight) {
     1529                if (variance) {
    15271530                    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);
    15311535                }
    15321536            } // Iterating over cell
  • branches/pap_branch_20090128/psModules/src/camera/pmFPACopy.c

    r20634 r21211  
    224224        targetReadout->data_exists = sourceReadout->data_exists;
    225225
    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);
    230231
    231232        // Copy bias
     
    328329        int xSize = psMetadataLookupS32(NULL, source->concepts, "CELL.XSIZE"); // CELL.XSIZE of source
    329330
    330         psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source");
     331        psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source");
    331332
    332333        if (sourceBin == 0) {
     
    478479        pmCell *targetCell = pmCellAlloc (targetChip, cellName);
    479480        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");
    481482        psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "", xParityTarget);
    482483        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");
    484485        psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "", yParityTarget);
    485486        if (!cellCopy(targetCell, sourceCell, true, 1, 1)) {
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAExtent.c

    r20635 r21211  
    1616    }
    1717    if (!image) {
    18         image = readout->weight;
     18        image = readout->variance;
    1919    }
    2020
     
    2323
    2424    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 {
    3939        // 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        }
    4343        return psRegionAlloc(0, xSize, 0, ySize);
    4444    }
     
    7070    // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent)
    7171    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;
    7878    }
    7979
     
    9797    // CELL.X0,Y0 are the coordinate of the amp on the chip, subtract size if parity flipped
    9898    if (xParityCell > 0) {
    99         cellExtent->x0 += x0;
    100         cellExtent->x1 += x0;
     99        cellExtent->x0 += x0;
     100        cellExtent->x1 += x0;
    101101    } 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;
    106106    }
    107107    if (yParityCell > 0) {
    108         cellExtent->y0 += y0;
    109         cellExtent->y1 += y0;
     108        cellExtent->y0 += y0;
     109        cellExtent->y1 += y0;
    110110    } 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;
    115115    }
    116116
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAMaskWeight.c

    r21183 r21211  
    4949}
    5050
    51 // Create the parent weight images that reside in the HDU
    52 static void createParentWeights(pmHDU *hdu // The HDU for which to create
     51// Create the parent variance images that reside in the HDU
     52static void createParentVariances(pmHDU *hdu // The HDU for which to create
    5353                               )
    5454{
     
    5858    // Generate the parent mask images
    5959    psArray *images = hdu->images;      // Array of images
    60     psArray *weights = hdu->weights;    // Array of weight images
    61     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;
    6464    }
    6565
    6666    for (long i = 0; i < images->n; i++) {
    6767        psImage *image = images->data[i]; // The image for this readout
    68         if (!image || weights->data[i]) {
     68        if (!image || variances->data[i]) {
    6969            continue;
    7070        }
    71         weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     71        variances->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    7272    }
    7373
     
    199199}
    200200
    201 bool pmReadoutSetWeight(pmReadout *readout, bool poisson)
     201bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
    202202{
    203203    PS_ASSERT_PTR_NON_NULL(readout, false);
     
    209209    float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
    210210    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");
    212212        return false;
    213213    }
    214214    float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
    215215    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");
    217217        return false;
    218218    }
     
    223223
    224224    if (poisson) {
    225         // Set weight image to the variance in ADU = f/g + rn^2
     225        // Set variance image to the variance in ADU = f/g + rn^2
    226226        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 weight must 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",
    232232                                               psScalarAlloc(1, PS_TYPE_F32));
    233233    } else {
    234234        // 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, "+",
    242242                                           psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
    243243
     
    245245}
    246246
    247 // this function creates the weight pixels, or uses the existing weight pixels.  it will set
    248 // the noise pixel values only if the weight image is not supplied
    249 bool pmReadoutGenerateWeight(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
     249bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson)
    250250{
    251251    PS_ASSERT_PTR_NON_NULL(readout, false);
     
    254254    bool mdok = true;                   // Status of MD lookup
    255255
    256     // Create the weight image if required
    257     if (readout->weight)
     256    // Create the variance image if required
     257    if (readout->variance)
    258258        return true;
    259259
    260260    psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
    261261    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");
    263263        return false;
    264264    }
     
    271271    }
    272272
    273     createParentWeights(hdu);
     273    createParentVariances(hdu);
    274274
    275275    // Need to identify which readout we're working with....
     
    280280    }
    281281
    282     psImage *weight = psImageSubset(hdu->weights->data[index], *trimsec); // The weight pixels
    283     if (!weight) {
     282    psImage *variance = psImageSubset(hdu->variances->data[index], *trimsec); // The variance pixels
     283    if (!variance) {
    284284        psString trimsecString = psRegionToString(*trimsec);
    285         psError(PS_ERR_UNKNOWN, false, "Unable to set weight from HDU with trimsec: %s.\n",
     285        psError(PS_ERR_UNKNOWN, false, "Unable to set variance from HDU with trimsec: %s.\n",
    286286                trimsecString);
    287287        psFree(trimsecString);
    288288        return false;
    289289    }
    290     psImageInit(weight, 0);
    291     readout->weight = weight;
    292 
    293     return pmReadoutSetWeight(readout, poisson);
    294 }
    295 
    296 bool pmReadoutGenerateMaskWeight(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     290    psImageInit(variance, 0);
     291    readout->variance = variance;
     292
     293    return pmReadoutSetVariance(readout, poisson);
     294}
     295
     296bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
    297297{
    298298    PS_ASSERT_PTR_NON_NULL(readout, false);
     
    301301
    302302    success &= pmReadoutGenerateMask(readout, satMask, badMask);
    303     success &= pmReadoutGenerateWeight(readout, poisson);
     303    success &= pmReadoutGenerateVariance(readout, poisson);
    304304
    305305    return success;
    306306}
    307307
    308 bool pmCellGenerateMaskWeight(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     308bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
    309309{
    310310    PS_ASSERT_PTR_NON_NULL(cell, false);
     
    314314    for (int i = 0; i < readouts->n; i++) {
    315315        pmReadout *readout = readouts->data[i]; // The readout
    316         success &= pmReadoutGenerateMaskWeight(readout, poisson, satMask, badMask);
     316        success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask);
    317317    }
    318318
     
    321321
    322322
    323 bool pmReadoutWeightRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
     323bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
    324324                                 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
    325325{
    326326    PM_ASSERT_READOUT_NON_NULL(readout, false);
    327327    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 parts
     328    PM_ASSERT_READOUT_VARIANCE(readout, false);
     329
     330    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
    331331
    332332    if (!psMemIncrRefCounter(rng)) {
     
    334334    }
    335335
    336     psStats *weightStats = psStatsAlloc(meanStat);// Statistics for mean
    337     if (!psImageBackground(weightStats, NULL, weight, mask, maskVal, rng)) {
    338         psError(PS_ERR_UNKNOWN, false, "Unable to measure mean weight for 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);
    340340        psFree(rng);
    341341        return false;
    342342    }
    343     float meanVariance = weightStats->robustMedian; // Mean variance
    344     psFree(weightStats);
     343    float meanVariance = varianceStats->robustMedian; // Mean variance
     344    psFree(varianceStats);
    345345
    346346    psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
     
    356356
    357357    float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be
    358     psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map 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));
    360360
    361361    return true;
     
    363363
    364364
    365 bool pmReadoutWeightRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
     365bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
    366366                               psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
    367367{
    368368    PM_ASSERT_READOUT_NON_NULL(readout, false);
    369369    PM_ASSERT_READOUT_IMAGE(readout, false);
    370     PM_ASSERT_READOUT_WEIGHT(readout, false);
     370    PM_ASSERT_READOUT_VARIANCE(readout, false);
    371371
    372372    if (!psMemIncrRefCounter(rng)) {
     
    374374    }
    375375
    376     psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
     376    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
    377377
    378378    // Measure background
     
    424424        float sumNoise = 0.0;       // Sum for noise measurement
    425425        float sumSource = 0.0;      // Sum for source measurement
    426         float sumWeight = 0.0;      // Sum for weight measurement
     426        float sumVariance = 0.0;      // Sum for variance measurement
    427427        float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels
    428428        for (int v = 0, y = yPix - size; v < fullSize; v++, y++) {
    429429            float xSumNoise = 0.0;  // Sum for noise measurement in x
    430430            float xSumSource = 0.0; // Sum for source measurement in x
    431             float xSumWeight = 0.0; // Sum for weight measurement in x
     431            float xSumVariance = 0.0; // Sum for variance measurement in x
    432432            float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x
    433433            float yGauss = gauss->data.F32[v]; // Value of Gaussian in y
     
    441441                xSumNoise += value * xGauss;
    442442                xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss;
    443                 xSumWeight += weight->data.F32[y][x] * xGauss2;
     443                xSumVariance += variance->data.F32[y][x] * xGauss2;
    444444                xSumGauss += xGauss;
    445445                xSumGauss2 += xGauss2;
     
    448448            sumNoise += xSumNoise * yGauss;
    449449            sumSource += xSumSource * yGauss;
    450             sumWeight += xSumWeight * yGauss2;
     450            sumVariance += xSumVariance * yGauss2;
    451451            sumGauss += xSumGauss * yGauss;
    452452            sumGauss2 += xSumGauss2 * yGauss2;
     
    454454
    455455        photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
    456                                                 isfinite(sumWeight) && sumGauss > 0 && sumGauss2 > 0) ?
     456                                                isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ?
    457457                                               0 : 0xFF);
    458458
    459459        float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise
    460460        float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source
    461         float smoothWeight = sumWeight / sumGauss2; // Value of smoothed weight pixel
     461        float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel
    462462
    463463        noise->data.F32[i] = smoothImageNoise;
    464464        source->data.F32[i] = smoothImageSource;
    465         guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothWeight : 0.0;
     465        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0;
    466466        psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
    467                 i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothWeight);
     467                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
    468468    }
    469469    psFree(gauss);
     
    516516    psFree(photMask);
    517517
    518     psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map 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));
    520520
    521521    return true;
     
    523523
    524524
    525 bool pmReadoutWeightRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
     525bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
    526526                           psStatsOptions stdevStat, int width, psRandom *rng)
    527527{
    528528    PM_ASSERT_READOUT_NON_NULL(readout, false);
    529529    PM_ASSERT_READOUT_IMAGE(readout, false);
    530     PM_ASSERT_READOUT_WEIGHT(readout, false);
     530    PM_ASSERT_READOUT_VARIANCE(readout, false);
    531531    PS_ASSERT_INT_POSITIVE(width, false);
    532532
     
    535535    }
    536536
    537     psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
     537    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
    538538    int numCols = image->numCols, numRows = image->numRows; // Size of images
    539539    int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions
     
    554554            psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
    555555            psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels
    556             psImage *subWeight = psImageSubset(weight, region); // Sub image of the weight pixels
     556            psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels
    557557            psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels
    558558
    559559            if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) ||
    560                 !psImageBackground(meanStats, &buffer, subWeight, subMask, maskVal, rng)) {
     560                !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) {
    561561                // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it.
    562562                psString regionStr = psRegionToString(region); // String with region
     
    564564                psFree(regionStr);
    565565                psErrorClear();
    566                 psImageInit(subWeight, NAN);
     566                psImageInit(subVariance, NAN);
    567567                if (subMask) {
    568568                    psImageInit(subMask, maskVal);
     
    574574                        "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n",
    575575                        xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar);
    576                 psBinaryOp(subWeight, subWeight, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
     576                psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
    577577            }
    578578
    579579            psFree(subImage);
    580             psFree(subWeight);
     580            psFree(subVariance);
    581581            psFree(subMask);
    582582        }
     
    597597
    598598    psImage *image = readout->image;    // Readout's image
    599     psImage *weight = readout->weight;  // Readout's weight
     599    psImage *variance = readout->variance;  // Readout's variance
    600600    int numCols = image->numCols, numRows = image->numRows; // Size of image
    601601
     
    607607    for (int y = 0; y < numRows; y++) {
    608608        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]))) {
    610610                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
    611611            }
     
    627627    psImageMaskType **maskData = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Dereference mask
    628628    psF32 **imageData = readout->image->data.F32;// Dereference image
    629     psF32 **weightData = readout->weight ? readout->weight->data.F32 : NULL; // Dereference weight map
     629    psF32 **varianceData = readout->variance ? readout->variance->data.F32 : NULL; // Dereference variance map
    630630
    631631    for (int y = 0; y < numRows; y++) {
     
    633633            if (maskData[y][x] & maskVal) {
    634634                imageData[y][x] = NAN;
    635                 if (weightData) {
    636                     weightData[y][x] = NAN;
     635                if (varianceData) {
     636                    varianceData[y][x] = NAN;
    637637                }
    638638            }
     
    656656    psImage *image = readout->image;    // Image
    657657    psImage *mask = readout->mask;      // Mask
    658     psImage *weight = readout->weight;  // Weight map
    659 
    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,
    661661                                                             NAN, NAN, maskBad, maskPoor, poorFrac, 0);
    662662    interp->shifting = false;           // Turn off "exact shifts" so we get proper interpolation
     
    666666    psPixels *pixels = psPixelsAllocEmpty(PIXELS_BUFFER); // Pixels that have been interpolated
    667667    psVector *imagePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for image
    668     psVector *weightPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for weight
     668    psVector *variancePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for variance
    669669    psVector *maskPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_IMAGE_MASK); // Corresponding values for mask
    670670    // NOTE: maskPix carries the actual image mask values -- do NOT use
     
    675675        for (int x = 0; x < numCols; x++) {
    676676            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
    677                 double imageValue, weightValue; // Image and weight value from interpolation
     677                double imageValue, varianceValue; // Image and variance value from interpolation
    678678                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);
    680680                if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) {
    681681                    psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y);
     
    683683                    psFree(pixels);
    684684                    psFree(imagePix);
    685                     psFree(weightPix);
     685                    psFree(variancePix);
    686686                    psFree(maskPix);
    687687                    return false;
     
    694694                pixels = psPixelsAdd(pixels, PIXELS_BUFFER, x, y);
    695695                imagePix = psVectorExtend(imagePix, PIXELS_BUFFER, 1);
    696                 weightPix = psVectorExtend(weightPix, PIXELS_BUFFER, 1);
     696                variancePix = psVectorExtend(variancePix, PIXELS_BUFFER, 1);
    697697                maskPix = psVectorExtend(maskPix, PIXELS_BUFFER, 1);
    698698                imagePix->data.F32[numBad] = imageValue;
    699                 weightPix->data.F32[numBad] = weightValue;
     699                variancePix->data.F32[numBad] = varianceValue;
    700700                maskPix->data.PS_TYPE_IMAGE_MASK_DATA[numBad] = maskValue;
    701701                numBad++;
     
    709709        int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of pixel
    710710        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];
    712712        mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskPix->data.PS_TYPE_IMAGE_MASK_DATA[i];
    713713    }
     
    715715    psFree(pixels);
    716716    psFree(imagePix);
    717     psFree(weightPix);
     717    psFree(variancePix);
    718718    psFree(maskPix);
    719719
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAMaskWeight.h

    r21183 r21211  
    55 * @author Eugene Magnier, IfA
    66 *
    7  * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2009-01-27 06:39:38 $
     7 * @version $Revision: 1.17.2.1 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-01-29 00:33:51 $
    99 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    1010 */
     
    4747
    4848
    49 /// Set a temporary readout weight map using CELL.GAIN and CELL.READNOISE
     49/// Set a temporary readout variance map using CELL.GAIN and CELL.READNOISE
    5050///
    51 /// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain
    52 /// (CELL.GAIN) and read noise (CELL.READNOISE).  The weight map that is produced within the readout is
    53 /// temporary --- it is not added to the HDU.  This is intended for when the user is iterating using
    54 /// pmReadoutReadNext, in which case the HDU can't be generated.
    55 bool pmReadoutSetWeight(pmReadout *readout, ///< Readout for which to set weight
    56                         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.
     55bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance
     56                          bool poisson    ///< Include poisson variance (in addition to read noise)?
     57    );
    5858
    5959/// Generate a readout mask (suitable for output) using CELL.SATURATION and CELL.BAD
     
    6666    );
    6767
    68 /// Generate a weight map (suitable for output) using CELL.GAIN and CELL.READNOISE
     68/// Generate a variance map (suitable for output) using CELL.GAIN and CELL.READNOISE
    6969///
    70 /// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain
    71 /// (CELL.GAIN) and read noise (CELL.READNOISE).  The weight map that is produced within the readout is
    72 /// suitable for output (complete with HDU entry).  This is intended for most operations.
    73 bool pmReadoutGenerateWeight(pmReadout *readout, ///< Readout for which to generate weight
    74                              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.
     73bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance
     74                               bool poisson    ///< Include poisson variance (in addition to read noise)?
     75    );
    7676
    77 /// Generate mask and weight map for a readout
     77/// Generate mask and variance map for a readout
    7878///
    79 /// Calls pmReadoutGenerateMask and pmReadoutGenerateWeight for the readout
    80 bool pmReadoutGenerateMaskWeight(pmReadout *readout, ///< Readout for which to generate mask and weights
    81                                  psImageMaskType sat, ///< Mask value to give saturated pixels
    82                                  psImageMaskType bad, ///< Mask value to give bad (low) pixels
    83                                  bool poisson ///< Use poisson weights (in addition to read noise)?
    84                                 );
     79/// Calls pmReadoutGenerateMask and pmReadoutGenerateVariance for the readout
     80bool 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    );
    8585
    86 /// Generate mask and weight maps for all readouts within a cell
     86/// Generate mask and variance maps for all readouts within a cell
    8787///
    88 /// Calls pmReadoutGenerateMaskWeight for each readout within the cell.
    89 bool pmCellGenerateMaskWeight(pmCell *cell, ///< Cell for which to generate mask and weights
    90                               psImageMaskType sat, ///< Mask value to give saturated pixels
    91                               psImageMaskType bad, ///< Mask value to give bad (low) pixels
    92                               bool poisson ///< Use poisson weights (in addition to read noise)?
    93                              );
     88/// Calls pmReadoutGenerateMaskVariance for each readout within the cell.
     89bool 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    );
    9494
    95 /// Renormalise the weight map to match the actual pixel variance
     95/// Renormalise the variance map to match the actual pixel variance
    9696///
    97 /// The weight (variance) map is adjusted so that the mean matches the actual pixel variance in the image
    98 bool pmReadoutWeightRenormPixels(
     97/// The variance map is adjusted so that the mean matches the actual pixel variance in the image
     98bool pmReadoutVarianceRenormPixels(
    9999    const pmReadout *readout,           ///< Readout to normalise
    100100    psImageMaskType maskVal,                 ///< Value to mask
     
    104104    );
    105105
    106 /// Renormalise the weight map to match the actual photometry variance
     106/// Renormalise the variance map to match the actual photometry variance
    107107///
    108 /// The weight (variance) map is adjusted so that the actual significance of fake sources matches the
     108/// The variance map is adjusted so that the actual significance of fake sources matches the
    109109/// guestimated significance
    110 bool pmReadoutWeightRenormPhot(
     110bool pmReadoutVarianceRenormPhot(
    111111    const pmReadout *readout,           ///< Readout to normalise
    112     psImageMaskType maskVal,                 ///< Value to mask
     112    psImageMaskType maskVal,            ///< Value to mask
    113113    int num,                            ///< Number of instances to measure over the image
    114114    float width,                        ///< Photometry width
     
    118118    );
    119119
    120 /// Renormalise the weight map to match the actual variance
     120/// Renormalise the variance map to match the actual variance
    121121///
    122122/// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for
    123123/// that patch corresponds.
    124 bool pmReadoutWeightRenorm(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)
     124bool 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)
    130130    );
    131131
     
    133133///
    134134/// 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 weight have their mask OR-ed with
     135/// necessary to mask them explicitly.  Non-finite pixels in the image or variance have their mask OR-ed with
    136136/// the provided value.
    137137bool pmReadoutMaskNonfinite(pmReadout *readout, ///< Readout to mask
     
    139139    );
    140140
    141 /// Apply a mask to the image and weight map
     141/// Apply a mask to the image and variance map
    142142///
    143143/// Unfortunately, image subtraction may result in a bi-modal image in masked areas, which can upset image
    144144/// 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.
    146146bool pmReadoutMaskApply(pmReadout *readout, ///< Readout to mask
    147147                        psImageMaskType maskVal ///< Mask value for which to apply mask
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAMosaic.c

    r21183 r21211  
    554554        // We have to do all of the hard work ourselves
    555555        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);
    566566          default:
    567567            psAbort("Should never get here.\n");
     
    575575static bool addCell(psArray *images,    // Array of images
    576576                    psArray *masks,     // Array of masks
    577                     psArray *weights,   // Array of weights
     577                    psArray *variances,   // Array of variances
    578578                    psVector *x0,       // Array of X0
    579579                    psVector *y0,       // Array of Y0
     
    604604        images  = psArrayRealloc(images,  index + CELL_LIST_BUFFER);
    605605        masks   = psArrayRealloc(masks,   index + CELL_LIST_BUFFER);
    606         weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER);
     606        variances = psArrayRealloc(variances, index + CELL_LIST_BUFFER);
    607607        x0    = psVectorRealloc(x0,    index+ CELL_LIST_BUFFER);
    608608        y0    = psVectorRealloc(y0,    index+ CELL_LIST_BUFFER);
     
    615615    images->n = index + 1;
    616616    masks->n = index + 1;
    617     weights->n = index + 1;
     617    variances->n = index + 1;
    618618    x0->n = index + 1;
    619619    y0->n = index + 1;
     
    726726    // The images to put into the mosaic
    727727    images->data[index]  = psMemIncrRefCounter(readout->image);
    728     weights->data[index] = psMemIncrRefCounter(readout->weight);
     728    variances->data[index] = psMemIncrRefCounter(readout->variance);
    729729    masks->data[index]   = psMemIncrRefCounter(readout->mask);
    730730
     
    740740static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
    741741                       psImage **mosaicMask, // The mosaic mask, to be returned
    742                        psImage **mosaicWeights, // The mosaic weights, to be returned
     742                       psImage **mosaicVariances, // The mosaic variances, to be returned
    743743                       int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
    744744                       const pmChip *chip, // Chip to mosaic
     
    749749    assert(mosaicImage);
    750750    assert(mosaicMask);
    751     assert(mosaicWeights);
     751    assert(mosaicVariances);
    752752    assert(xBinChip);
    753753    assert(yBinChip);
     
    756756
    757757    psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
    758     psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
     758    psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked
    759759    psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
    760760    psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
     
    802802            continue;
    803803        }
    804         allGood &= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
     804        allGood &= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip,
    805805                           cell, xBinChip, yBinChip, false, x0Target, y0Target,
    806806                           xParityTarget, yParityTarget);
     
    826826    if (allGood) {
    827827        *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
    828         *mosaicWeights = 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);
    829829        *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank);
    830830    }
     
    832832    // Clean up
    833833    psFree(images);
    834     psFree(weights);
     834    psFree(variances);
    835835    psFree(masks);
    836836    psFree(xFlip);
     
    847847static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
    848848                      psImage **mosaicMask, // The mosaic mask, to be returned
    849                       psImage **mosaicWeights, // The mosaic weights, to be returned
     849                      psImage **mosaicVariances, // The mosaic variances, to be returned
    850850                      int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
    851851                      const pmFPA *fpa,  // FPA to mosaic
     
    857857    assert(mosaicImage);
    858858    assert(mosaicMask);
    859     assert(mosaicWeights);
     859    assert(mosaicVariances);
    860860    assert(xBinFPA);
    861861    assert(yBinFPA);
     
    865865
    866866    psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
    867     psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
     867    psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked
    868868    psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
    869869    psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
     
    941941                continue;
    942942            }
    943             allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
     943            allGood |= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip,
    944944                               cell, xBinFPA, yBinFPA, true, x0Target, y0Target,
    945945                               xParityTarget, yParityTarget);
     
    960960    if (allGood) {
    961961        *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
    962         *mosaicWeights = 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);
    963963        *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank);
    964964    }
     
    966966    // Clean up
    967967    psFree(images);
    968     psFree(weights);
     968    psFree(variances);
    969969    psFree(masks);
    970970    psFree(xFlip);
     
    10251025    psImage *mosaicImage   = NULL;      // The mosaic image
    10261026    psImage *mosaicMask    = NULL;      // The mosaic mask
    1027     psImage *mosaicWeights = NULL;      // The mosaic weights
     1027    psImage *mosaicVariances = NULL;      // The mosaic variances
    10281028
    10291029    // Find the HDU
     
    10511051            }
    10521052        }
    1053         if (hdu->weights) {
    1054             mosaicWeights = psImageSubset(hdu->weights->data[0], bounds);
    1055             if (!mosaicWeights) {
    1056                 psError(PS_ERR_UNKNOWN, false, "Unable to select weight pixels.\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");
    10571057                return false;
    10581058            }
     
    10611061        // Case 2 --- we need to mosaic by cut and paste
    10621062        psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
    1063         if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source, targetCell, blank)) {
     1063        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) {
    10641064            psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
    10651065            return false;
     
    10941094    newReadout->image  = mosaicImage;
    10951095    newReadout->mask   = mosaicMask;
    1096     newReadout->weight = mosaicWeights;
     1096    newReadout->variance = mosaicVariances;
    10971097    psFree(newReadout);                 // Drop reference
    10981098
     
    12621262    psImage *mosaicImage   = NULL;      // The mosaic image
    12631263    psImage *mosaicMask    = NULL;      // The mosaic mask
    1264     psImage *mosaicWeights = NULL;      // The mosaic weights
     1264    psImage *mosaicVariances = NULL;      // The mosaic variances
    12651265
    12661266    // Find the HDU
     
    12751275            mosaicMask = psImageSubset(hdu->masks->data[0], *fpaRegion);
    12761276        }
    1277         if (hdu->weights) {
    1278             mosaicWeights = psImageSubset(hdu->weights->data[0], *fpaRegion);
     1277        if (hdu->variances) {
     1278            mosaicVariances = psImageSubset(hdu->variances->data[0], *fpaRegion);
    12791279        }
    12801280    } else {
    12811281        // Case 2 --- we need to mosaic by cut and paste
    12821282        psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
    1283         if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source,
     1283        if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source,
    12841284                       targetChip, targetCell, blank)) {
    12851285            psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n");
     
    13381338    newReadout->image  = mosaicImage;
    13391339    newReadout->mask   = mosaicMask;
    1340     newReadout->weight = mosaicWeights;
     1340    newReadout->variance = mosaicVariances;
    13411341    psFree(newReadout);                 // Drop reference
    13421342
  • branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c

    r21183 r21211  
    2828    FPA_READ_TYPE_IMAGE,                // Read image
    2929    FPA_READ_TYPE_MASK,                 // Read mask
    30     FPA_READ_TYPE_WEIGHT,               // Read weight map
     30    FPA_READ_TYPE_VARIANCE,               // Read variance map
    3131    FPA_READ_TYPE_HEADER                // Read header
    3232} fpaReadType;
     
    5454      case FPA_READ_TYPE_MASK:
    5555        return readout->thisMaskScan;
    56       case FPA_READ_TYPE_WEIGHT:
    57         return readout->thisWeightScan;
     56      case FPA_READ_TYPE_VARIANCE:
     57        return readout->thisVarianceScan;
    5858      default:
    5959        psAbort("Unknown read type: %x\n", type);
     
    7474        readout->thisMaskScan = thisScan;
    7575        return readout->lastMaskScan;
    76       case FPA_READ_TYPE_WEIGHT:
    77         readout->thisWeightScan = thisScan;
    78         return readout->lastWeightScan;
     76      case FPA_READ_TYPE_VARIANCE:
     77        readout->thisVarianceScan = thisScan;
     78        return readout->lastVarianceScan;
    7979      default:
    8080        psAbort("Unknown read type: %x\n", type);
     
    9393      case FPA_READ_TYPE_MASK:
    9494        return readout->lastMaskScan;
    95       case FPA_READ_TYPE_WEIGHT:
    96         return readout->lastWeightScan;
     95      case FPA_READ_TYPE_VARIANCE:
     96        return readout->lastVarianceScan;
    9797      default:
    9898        psAbort("Unknown read type: %x\n", type);
     
    113113        readout->lastMaskScan = lastScan;
    114114        return readout->lastMaskScan;
    115       case FPA_READ_TYPE_WEIGHT:
    116         readout->lastWeightScan = lastScan;
    117         return readout->lastWeightScan;
     115      case FPA_READ_TYPE_VARIANCE:
     116        readout->lastVarianceScan = lastScan;
     117        return readout->lastVarianceScan;
    118118      default:
    119119        psAbort("Unknown read type: %x\n", type);
     
    132132      case FPA_READ_TYPE_MASK:
    133133        return &readout->mask;
    134       case FPA_READ_TYPE_WEIGHT:
    135         return &readout->weight;
     134      case FPA_READ_TYPE_VARIANCE:
     135        return &readout->variance;
    136136      default:
    137137        psAbort("Unknown read type: %x\n", type);
     
    350350    *target = psImageSubset(image, region);
    351351
    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)
    353353    if (type == FPA_READ_TYPE_IMAGE) {
    354354        if (readout->bias->n != 0) {
     
    582582    *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]);
    583583
    584     // Read overscans only for "image" type --- weights and masks shouldn't record overscans
     584    // Read overscans only for "image" type --- variances and masks shouldn't record overscans
    585585    if (type == FPA_READ_TYPE_IMAGE) {
    586586        // Blow away existing data
     
    608608}
    609609
    610 // Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadWeight
     610// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadVariance
    611611// Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts.
    612612static bool cellRead(pmCell *cell,      // Cell into which to read
     
    640640        dataPointer = hdu->masks;
    641641        break;
    642       case FPA_READ_TYPE_WEIGHT:
    643         hduReadFunc = pmHDUReadWeight;
    644         dataPointer = hdu->weights;
     642      case FPA_READ_TYPE_VARIANCE:
     643        hduReadFunc = pmHDUReadVariance;
     644        dataPointer = hdu->variances;
    645645        break;
    646646      default:
     
    680680        imageArray = hdu->masks;
    681681        break;
    682       case FPA_READ_TYPE_WEIGHT:
    683         imageArray = hdu->weights;
     682      case FPA_READ_TYPE_VARIANCE:
     683        imageArray = hdu->variances;
    684684        break;
    685685      default:
     
    729729
    730730
    731 // Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadWeight
     731// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadVariance
    732732// Iterates over component cells, reading each
    733733static bool chipRead(pmChip *chip,      // Chip into which to read
     
    760760
    761761
    762 // Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadWeight
     762// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadVariance
    763763// Iterates over component chips, reading each
    764764static bool fpaRead(pmFPA *fpa,         // FPA into which to read
     
    10911091
    10921092//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1093 // Reading the weight map
    1094 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1095 
    1096 bool pmReadoutMoreWeight(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
     1093// Reading the variance map
     1094//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     1095
     1096bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
    10971097{
    10981098    PS_ASSERT_PTR_NON_NULL(readout, false);
    10991099    PS_ASSERT_FITS_NON_NULL(fits, false);
    11001100
    1101     return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_WEIGHT, config);
    1102 }
    1103 
    1104 bool pmReadoutReadChunkWeight(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
     1101    return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config);
     1102}
     1103
     1104bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
    11051105                              pmConfig *config)
    11061106{
     
    11101110    PS_ASSERT_INT_NONNEGATIVE(numScans, false);
    11111111
    1112     return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_WEIGHT, config);
    1113 }
    1114 
    1115 bool pmReadoutReadWeight(pmReadout *readout, psFits *fits, int z, pmConfig *config)
     1112    return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config);
     1113}
     1114
     1115bool pmReadoutReadVariance(pmReadout *readout, psFits *fits, int z, pmConfig *config)
    11161116{
    11171117    PS_ASSERT_PTR_NON_NULL(readout, false);
    11181118    PS_ASSERT_FITS_NON_NULL(fits, false);
    11191119
    1120     return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_WEIGHT, config);
    1121 }
    1122 
    1123 bool pmCellReadWeight(pmCell *cell, psFits *fits, pmConfig *config)
     1120    return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config);
     1121}
     1122
     1123bool pmCellReadVariance(pmCell *cell, psFits *fits, pmConfig *config)
    11241124{
    11251125    PS_ASSERT_PTR_NON_NULL(cell, false);
    11261126    PS_ASSERT_FITS_NON_NULL(fits, false);
    11271127
    1128     return cellRead(cell, fits, config, FPA_READ_TYPE_WEIGHT);
    1129 }
    1130 
    1131 bool pmChipReadWeight(pmChip *chip, psFits *fits, pmConfig *config)
     1128    return cellRead(cell, fits, config, FPA_READ_TYPE_VARIANCE);
     1129}
     1130
     1131bool pmChipReadVariance(pmChip *chip, psFits *fits, pmConfig *config)
    11321132{
    11331133    PS_ASSERT_PTR_NON_NULL(chip, false);
    11341134    PS_ASSERT_FITS_NON_NULL(fits, false);
    11351135
    1136     return chipRead(chip, fits, config, FPA_READ_TYPE_WEIGHT);
    1137 }
    1138 
    1139 bool pmFPAReadWeight(pmFPA *fpa, psFits *fits, pmConfig *config)
     1136    return chipRead(chip, fits, config, FPA_READ_TYPE_VARIANCE);
     1137}
     1138
     1139bool pmFPAReadVariance(pmFPA *fpa, psFits *fits, pmConfig *config)
    11401140{
    11411141    PS_ASSERT_PTR_NON_NULL(fpa, false);
    11421142    PS_ASSERT_FITS_NON_NULL(fits, false);
    11431143
    1144     return fpaRead(fpa, fits, config, FPA_READ_TYPE_WEIGHT);
     1144    return fpaRead(fpa, fits, config, FPA_READ_TYPE_VARIANCE);
    11451145}
    11461146
  • branches/pap_branch_20090128/psModules/src/camera/pmFPARead.h

    r18163 r21211  
    44 * @author Paul Price, IfA
    55 *
    6  * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2008-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 $
    88 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    99 */
     
    6060                       int numRows,      // The number of rows to read
    6161                       pmConfig *config
    62                       );
     62    );
    6363
    6464/// Return the number of readouts within a cell
    6565///
    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).
    6767int pmCellNumReadouts(pmCell *cell, psFits *fits, pmConfig *config);
    6868
     
    7474                psFits *fits,           // FITS file from which to read
    7575                pmConfig *config        // Configuration
    76                );
     76    );
    7777
    7878/// Read a chip
     
    8282                psFits *fits,           // FITS file from which to read
    8383                pmConfig *config        // Configuration
    84                );
     84    );
    8585
    8686/// Read an FPA
     
    9090               psFits *fits,            // FITS file from which to read
    9191               pmConfig *config         // Configuration
    92               );
     92    );
    9393
    9494// Mask functions follow
     
    126126                    psFits *fits,       // FITS file from which to read
    127127                    pmConfig *config    // Configuration
    128                    );
     128    );
    129129
    130130/// Read an entire chip into the mask
     
    134134                    psFits *fits,       // FITS file from which to read
    135135                    pmConfig *config    // Configuration
    136                    );
     136    );
    137137
    138138/// Read an entire FPA into the mask
     
    142142                   psFits *fits,        // FITS file from which to read
    143143                   pmConfig *config     // Configuration
    144                   );
    145 
    146 // Weight functions follow
    147 
    148 /// Check to see if there is more to read when reading a readout incrementally into the weight
    149 bool pmReadoutMoreWeight(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 weight
     144    );
     145
     146// Variance functions follow
     147
     148/// Check to see if there is more to read when reading a readout incrementally into the variance
     149bool 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
    157157///
    158158/// Allows reading the readout incrementally
    159 bool pmReadoutReadChunkWeight(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 weight
    168 bool pmReadoutReadWeight(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 weight
    175 ///
    176 /// Same as pmCellRead, but reads into the weight element of the readouts.
    177 bool pmCellReadWeight(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 weight
    183 ///
    184 /// Same as pmChipRead, but reads into the weight element of the readouts.
    185 bool pmChipReadWeight(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 weight
    191 ///
    192 /// Same as pmFPARead, but reads into the weight element of the readouts.
    193 bool pmFPAReadWeight(pmFPA *fpa,        // FPA to read into
    194                      psFits *fits,      // FITS file from which to read
    195                      pmConfig *config   // Configuration
    196                     );
     159bool 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
     168bool 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.
     177bool 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.
     185bool 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.
     193bool pmFPAReadVariance(pmFPA *fpa,        // FPA to read into
     194                       psFits *fits,      // FITS file from which to read
     195                       pmConfig *config   // Configuration
     196    );
    197197
    198198/// Read cell headers
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c

    r19385 r21211  
    2929    FPA_WRITE_TYPE_IMAGE,               // Write image
    3030    FPA_WRITE_TYPE_MASK,                // Write mask
    31     FPA_WRITE_TYPE_WEIGHT               // Write weight map
     31    FPA_WRITE_TYPE_VARIANCE               // Write variance map
    3232} fpaWriteType;
    3333
     
    4646    case FPA_WRITE_TYPE_MASK:
    4747        return &hdu->masks;
    48     case FPA_WRITE_TYPE_WEIGHT:
    49         return &hdu->weights;
     48    case FPA_WRITE_TYPE_VARIANCE:
     49        return &hdu->variances;
    5050    default:
    5151        psAbort("Unknown write type: %x\n", type);
     
    6666    case FPA_WRITE_TYPE_MASK:
    6767        return pmHDUWriteMask(hdu, fits, config);
    68     case FPA_WRITE_TYPE_WEIGHT:
    69         return pmHDUWriteWeight(hdu, fits, config);
     68    case FPA_WRITE_TYPE_VARIANCE:
     69        return pmHDUWriteVariance(hdu, fits, config);
    7070    default:
    7171        psAbort("Unknown write type: %x\n", type);
     
    7474}
    7575
    76 // Write a cell image/mask/weight
     76// Write a cell image/mask/variance
    7777static bool cellWrite(pmCell *cell,     // Cell to write
    7878                      psFits *fits,     // FITS file to which to write
     
    124124}
    125125
    126 // Write a chip image/mask/weight
     126// Write a chip image/mask/variance
    127127static bool chipWrite(pmChip *chip,     // Chip to write
    128128                      psFits *fits,     // FITS file to which to write
     
    188188
    189189
    190 // Write an FPA image/mask/weight
     190// Write an FPA image/mask/variance
    191191static bool fpaWrite(pmFPA *fpa,        // FPA to write
    192192                     psFits *fits,      // FITS file to which to write
     
    417417
    418418
    419 bool pmCellWriteWeight(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
     419bool pmCellWriteVariance(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
    420420{
    421421    PS_ASSERT_PTR_NON_NULL(cell, false);
    422422    PS_ASSERT_PTR_NON_NULL(fits, false);
    423     return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_WEIGHT);
    424 }
    425 
    426 bool pmChipWriteWeight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
     423    return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_VARIANCE);
     424}
     425
     426bool pmChipWriteVariance(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
    427427{
    428428    PS_ASSERT_PTR_NON_NULL(chip, false);
    429429    PS_ASSERT_PTR_NON_NULL(fits, false);
    430     return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT);
    431 }
    432 
    433 bool pmFPAWriteWeight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
     430    return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE);
     431}
     432
     433bool pmFPAWriteVariance(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
    434434{
    435435    PS_ASSERT_PTR_NON_NULL(fpa, false);
    436436    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);
    438438}
    439439
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.h

    r18163 r21211  
    44 * @author Paul Price, IfA
    55 *
    6  * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2008-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 $
    88 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    99 */
     
    2424                        psFits *fits,   ///<  FITS file to which to write
    2525                        int z           ///<  Image plane to write
    26                        );
     26    );
    2727
    2828/// Write a cell to a FITS file
     
    3636                 pmConfig *config,      ///<  Configuration
    3737                 bool blank             ///<  Write a blank PHU?
    38                 );
     38    );
    3939
    4040/// Write a chip to a FITS file
     
    4949                 bool blank,            ///<  Write a blank PHU?
    5050                 bool recurse           ///<  Recurse to lower levels?
    51                 );
     51    );
    5252
    5353/// Write an FPA to a FITS file
     
    6262                bool blank,             ///<  Write a blank PHU?
    6363                bool recurse            ///<  Recurse to lower levels?
    64                );
     64    );
    6565
    6666/// Write a cell mask to a FITS file
     
    7474                     pmConfig *config,  ///<  Configuration
    7575                     bool blank         ///<  Write a blank PHU?
    76                     );
     76    );
    7777
    7878/// Write a chip mask to a FITS file
     
    8888                     bool blank,        ///<  Write a blank PHU?
    8989                     bool recurse       ///<  Recurse to lower levels?
    90                     );
     90    );
    9191
    9292/// Write an FPA mask to a FITS file
     
    102102                    bool blank,         ///<  Write a blank PHU?
    103103                    bool recurse        ///<  Recurse to lower levels?
    104                    );
     104    );
    105105
    106 /// Write a cell weight to a FITS file
     106/// Write a cell variance to a FITS file
    107107///
    108 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     108/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    109109/// 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 loop
    111 /// with blank=true in order to produce the correct file structure.
    112 bool pmCellWriteWeight(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                       );
     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.
     112bool 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    );
    117117
    118 /// Write a chip weight to a FITS file
     118/// Write a chip variance to a FITS file
    119119///
    120 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     120/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    121121/// 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 at
    123 /// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct
     122/// 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
    124124/// file structure.
    125 bool pmChipWriteWeight(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                       );
     125bool 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    );
    131131
    132 /// Write an FPA weight to a FITS file
     132/// Write an FPA variance to a FITS file
    133133///
    134 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     134/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    135135/// 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 at
    137 /// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct
     136/// 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
    138138/// file structure.
    139 bool pmFPAWriteWeight(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                      );
     139bool 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    );
    145145
    146146
     
    153153                     const pmCell *cell, ///< Cell containing FITS table in the analysis metadata
    154154                     const char *name   ///< Name for the table data, and the extension name
    155                     );
     155    );
    156156
    157157int pmChipWriteTable(psFits *fits,      ///< FITS file to which to write
    158158                     const pmChip *chip, ///< Chip containing cells with tables to write
    159159                     const char *name   ///< Name for the table data, and the extension name
    160                     );
     160    );
    161161
    162162int pmFPAWriteTable(psFits *fits,       ///< FITS file to which to write
    163163                    const pmFPA *fpa,   ///< FPA containing cells with tables to write
    164164                    const char *name    ///< Name for the table data, and the extension name
    165                    );
     165    );
    166166
    167167// XXX better name, please
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAfile.c

    r19307 r21211  
    473473        return PM_FPA_FILE_MASK;
    474474    }
    475     if (!strcasecmp(type, "WEIGHT"))     {
    476         return PM_FPA_FILE_WEIGHT;
     475    if (!strcasecmp(type, "VARIANCE"))     {
     476        return PM_FPA_FILE_VARIANCE;
    477477    }
    478478    if (!strcasecmp(type, "FRINGE")) {
     
    527527      case PM_FPA_FILE_MASK:
    528528        return ("MASK");
    529       case PM_FPA_FILE_WEIGHT:
    530         return ("WEIGHT");
     529      case PM_FPA_FILE_VARIANCE:
     530        return ("VARIANCE");
    531531      case PM_FPA_FILE_FRINGE:
    532532        return ("FRINGE");
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAfile.h

    r19307 r21211  
    44 * @author EAM, IfA
    55 *
    6  * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2008-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 $
    88 * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
    99 */
     
    3737    PM_FPA_FILE_IMAGE,
    3838    PM_FPA_FILE_MASK,
    39     PM_FPA_FILE_WEIGHT,
     39    PM_FPA_FILE_VARIANCE,
    4040    PM_FPA_FILE_FRINGE,
    4141    PM_FPA_FILE_DARK,
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAfileFitsIO.c

    r18601 r21211  
    148148      case PM_FPA_FILE_IMAGE:
    149149      case PM_FPA_FILE_MASK:
    150       case PM_FPA_FILE_WEIGHT:
     150      case PM_FPA_FILE_VARIANCE:
    151151      case PM_FPA_FILE_HEADER:
    152152      case PM_FPA_FILE_FRINGE:
     
    329329}
    330330
    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);
     331bool 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);
    336337}
    337338
     
    347348    PS_ASSERT_PTR_NON_NULL(view, false);
    348349    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);
    350352}
    351353
     
    432434}
    433435
    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);
     436bool 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);
    439442}
    440443
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAfileFitsIO.h

    r18601 r21211  
    55 * @author PAP, IfA
    66 *
    7  * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2008-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 $
    99 * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
    1010 */
     
    2727                            pmConfig *config
    2828                          );
    29 /// Read a weight map into the current view
    30 bool pmFPAviewReadFitsWeight(const pmFPAview *view,  ///< View specifying level of interest
     29/// Read a variance map into the current view
     30bool pmFPAviewReadFitsVariance(const pmFPAview *view,  ///< View specifying level of interest
    3131                             pmFPAfile *file, ///< FPA file into which to read
    3232                            pmConfig *config
     
    5757                           );
    5858
    59 /// Write the weight map for the specified view
    60 bool pmFPAviewWriteFitsWeight(const pmFPAview *view, ///< View specifying level of interest
     59/// Write the variance map for the specified view
     60bool pmFPAviewWriteFitsVariance(const pmFPAview *view, ///< View specifying level of interest
    6161                              pmFPAfile *file, ///< FPA file to write
    6262                              pmConfig *config ///< Configuration
  • branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c

    r20937 r21211  
    173173        status = pmFPAviewReadFitsMask(view, file, config);
    174174        break;
    175       case PM_FPA_FILE_WEIGHT:
    176         status = pmFPAviewReadFitsWeight(view, file, config);
     175      case PM_FPA_FILE_VARIANCE:
     176        status = pmFPAviewReadFitsVariance(view, file, config);
    177177        break;
    178178      case PM_FPA_FILE_HEADER:
     
    262262      case PM_FPA_FILE_IMAGE:
    263263      case PM_FPA_FILE_MASK:
    264       case PM_FPA_FILE_WEIGHT:
     264      case PM_FPA_FILE_VARIANCE:
    265265      case PM_FPA_FILE_FRINGE:
    266266      case PM_FPA_FILE_DARK: {
     
    427427        status = pmFPAviewWriteFitsMask(view, file, config);
    428428        break;
    429       case PM_FPA_FILE_WEIGHT:
    430         status = pmFPAviewWriteFitsWeight(view, file, config);
     429      case PM_FPA_FILE_VARIANCE:
     430        status = pmFPAviewWriteFitsVariance(view, file, config);
    431431        break;
    432432      case PM_FPA_FILE_HEADER:
     
    524524      case PM_FPA_FILE_IMAGE:
    525525      case PM_FPA_FILE_MASK:
    526       case PM_FPA_FILE_WEIGHT:
     526      case PM_FPA_FILE_VARIANCE:
    527527      case PM_FPA_FILE_HEADER:
    528528      case PM_FPA_FILE_FRINGE:
     
    589589      case PM_FPA_FILE_IMAGE:
    590590      case PM_FPA_FILE_MASK:
    591       case PM_FPA_FILE_WEIGHT:
     591      case PM_FPA_FILE_VARIANCE:
    592592      case PM_FPA_FILE_HEADER:
    593593      case PM_FPA_FILE_FRINGE:
     
    730730    psString tmpName = pmConfigConvertFilename (file->filename, config, create, false);
    731731    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;
    734734    }
    735735    psFree (file->filename);
     
    740740      case PM_FPA_FILE_IMAGE:
    741741      case PM_FPA_FILE_MASK:
    742       case PM_FPA_FILE_WEIGHT:
     742      case PM_FPA_FILE_VARIANCE:
    743743      case PM_FPA_FILE_HEADER:
    744744      case PM_FPA_FILE_FRINGE:
     
    783783          // determine camera if not specified already
    784784          // 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;
    788788          file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true);
    789789          if (!file->format) {
     
    794794          psFree(phu);
    795795
    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;
    810810        }
    811811
     
    938938      case PM_FPA_FILE_IMAGE:
    939939      case PM_FPA_FILE_MASK:
    940       case PM_FPA_FILE_WEIGHT:
     940      case PM_FPA_FILE_VARIANCE:
    941941      case PM_FPA_FILE_DARK:
    942942      case PM_FPA_FILE_FRINGE:
  • branches/pap_branch_20090128/psModules/src/camera/pmHDU.c

    r21183 r21211  
    4545    psFree(hdu->header);
    4646    psFree(hdu->images);
    47     psFree(hdu->weights);
     47    psFree(hdu->variances);
    4848    psFree(hdu->masks);
    4949}
     
    6969    hdu->header  = NULL;
    7070    hdu->images  = NULL;
    71     hdu->weights = NULL;
     71    hdu->variances = NULL;
    7272    hdu->masks   = NULL;
    7373
     
    155155}
    156156
    157 bool pmHDUReadWeight(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);
     157bool 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);
    163163}
    164164
     
    235235}
    236236
    237 bool pmHDUWriteWeight(pmHDU *hdu, psFits *fits, const pmConfig *config)
     237bool pmHDUWriteVariance(pmHDU *hdu, psFits *fits, const pmConfig *config)
    238238{
    239239    PS_ASSERT_PTR_NON_NULL(hdu, false);
     
    241241
    242242    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  
    44 * @author Paul Price, IfA
    55 *
    6  * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2008-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 $
    88 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    99 */
     
    2020/// An instance of the FITS Header Data Unit
    2121///
    22 /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weight data, but these are
     22/// Of course, it is not an exact replica of a FITS HDU --- they have no mask and variance data, but these are
    2323/// stored here for convenience --- it keeps all the relevant data about the image in one place.
    2424typedef struct
     
    2929    psMetadata *header;                 ///< The FITS header, or NULL if primary for FITS; or section info
    3030    psArray *images;                    ///< The pixel data
    31     psArray *weights;                   ///< The pixel data
     31    psArray *variances;                   ///< The pixel data
    3232    psArray *masks;                     ///< The pixel data
    3333}
     
    6060                  );
    6161
    62 /// Read the HDU header and weight map
     62/// Read the HDU header and variance map
    6363///
    6464/// Moves to the appropriate extension
    65 bool pmHDUReadWeight(pmHDU *hdu,        ///< HDU to read
    66                      psFits *fits       ///< FITS file to read from
     65bool pmHDUReadVariance(pmHDU *hdu,        ///< HDU to read
     66                       psFits *fits       ///< FITS file to read from
    6767    );
    6868
     
    7979    );
    8080
    81 /// Write the HDU header and weight map
    82 bool pmHDUWriteWeight(pmHDU *hdu,       ///< HDU to write
    83                       psFits *fits,     ///< FITS file to write to
    84                       const pmConfig *config  ///< Configuration
     81/// Write the HDU header and variance map
     82bool pmHDUWriteVariance(pmHDU *hdu,       ///< HDU to write
     83                        psFits *fits,     ///< FITS file to write to
     84                        const pmConfig *config  ///< Configuration
    8585    );
    8686
  • branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c

    r18227 r21211  
    274274        pmReadout *readout = cell->readouts->data[0]; // The first readout, as representative
    275275        // 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);
    277277        if (!image) {
    278278            continue;
     
    283283                     image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows);
    284284        }
    285         if (readout->weight &&
    286                 (readout->weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) {
    287             psWarning("Image and weight have 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);
    289289        }
    290290        // New reference
     
    355355    psElemType imageType = 0;           // Type of readout images
    356356    psElemType maskType = 0;            // Type of readout masks
    357     psElemType weightType = 0;          // Type of readout weights
     357    psElemType varianceType = 0;          // Type of readout variances
    358358    {
    359359        psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells
     
    381381                    maskType = checkTypes(maskType, readout->mask->type.type);
    382382                }
    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);
    385385                }
    386386            }
     
    388388        psFree(iter);
    389389    }
    390     if (numReadouts == 0 || (imageType == 0 && maskType == 0 && weightType == 0)) {
     390    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) {
    391391        // Nothing from which to create an HDU
    392392        psFree(cells);
     
    421421        }
    422422    }
    423     if (weightType) {
    424         hdu->weights = psArrayAlloc(numReadouts);
     423    if (varianceType) {
     424        hdu->variances = psArrayAlloc(numReadouts);
    425425        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;
    429429        }
    430430    }
     
    450450            psArray *hduImages = hdu->images; // Array of images in the HDU
    451451            psArray *hduMasks = hdu->masks; // Array of masks in the HDU
    452             psArray *hduWeights = hdu->weights; // Array of weights in the HDU
     452            psArray *hduVariances = hdu->variances; // Array of variances in the HDU
    453453            for (int i = 0; i < readouts->n; i++) {
    454454                pmReadout *readout = readouts->data[i]; // The readout of interest
     
    467467                    readout->mask = new;
    468468                }
    469                 if (readout->weight) {
    470                     psImage *new = pasteImage(hduWeights->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;
    473473                }
    474474
     
    581581        return false;
    582582    }
    583     if (hdu->images && hdu->masks && hdu->weights) {
     583    if (hdu->images && hdu->masks && hdu->variances) {
    584584        // It's already here!
    585585        return true;
     
    631631        return generateForCells(chip);
    632632    }
    633     if (hdu->images && hdu->masks && hdu->weights) {
     633    if (hdu->images && hdu->masks && hdu->variances) {
    634634        // It's already here!
    635635        return true;
     
    680680        return generateForChips(fpa);
    681681    }
    682     if (hdu->images && hdu->masks && hdu->weights) {
     682    if (hdu->images && hdu->masks && hdu->variances) {
    683683        // It's already here!
    684684        return true;
  • branches/pap_branch_20090128/psModules/src/camera/pmHDUUtils.c

    r18554 r21211  
    1818
    1919    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        }
    2929    }
    3030    return NULL;
     
    158158
    159159    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]; // Weight image of interest
     160    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
    164164            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);
    166166        }
    167167    } else {
    168         fprintf(fd, "NO weights.\n");
     168        fprintf(fd, "NO variances.\n");
    169169    }
    170170
  • branches/pap_branch_20090128/psModules/src/camera/pmReadoutStack.c

    r21183 r21211  
    1212// XXX should it be an error for the image to exist?
    1313psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image
    14                                    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
     14                                   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
    1818    )
    1919{
     
    2424
    2525    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");
    2727    }
    2828    psImageInit(image, init);
     
    3535// XXX not sure why this should call psMemIncrRefCounter
    3636psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image
    37                                    const char *name       // Name of image in analysis metadata
     37                                   const char *name       // Name of image in analysis metadata
    3838    )
    3939{
     
    7878
    7979// 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)
     80bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool variance, psImageMaskType blank)
    8181{
    8282    PS_ASSERT_PTR_NON_NULL(readout, false);
     
    9292
    9393    if (mask) {
    94         // XXX is this an error?
     94        // XXX is this an error?
    9595        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);
    105105    }
    106106
     
    162162
    163163bool 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,
    165165                         psImageMaskType blank)
    166166{
     
    203203    }
    204204
    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 *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    212             psImageInit(newWeight, NAN);
    213             psImageOverlaySection(newWeight, 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;
    216216        }
    217217    }
  • branches/pap_branch_20090128/psModules/src/detrend/pmShutterCorrection.c

    r21183 r21211  
    383383    psArray *images = psArrayAlloc(num);// Array of images
    384384    psArray *masks = NULL; // Array of masks
    385     psArray *weights = NULL; // Array of weights
     385    psArray *variances = NULL; // Array of variances
    386386    psVector *exptimes = psVectorAlloc(num, PS_TYPE_F32); // Vector of exposure times
    387387
     
    392392            masks = psArrayAlloc(num);
    393393        }
    394         if (readout->weight)
     394        if (readout->variance)
    395395        {
    396             weights = psArrayAlloc(num);
     396            variances = psArrayAlloc(num);
    397397        }
    398398    }
     
    472472        }
    473473
    474         psImage *weight = readout->weight; // Weight map of interest
    475         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");
    478478                goto MEASURE_ERROR;
    479479            }
    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);
    484484                goto MEASURE_ERROR;
    485485            }
    486             if (weight->numRows != numRows || weight->numCols != numCols) {
     486            if (variance->numRows != numRows || variance->numCols != numCols) {
    487487                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,
    489489                        numCols, numRows);
    490490                goto MEASURE_ERROR;
    491491            }
    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");
    494494            goto MEASURE_ERROR;
    495495        }
     
    566566    psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef);
    567567
    568     // Check the weights
    569     if (weights && nIter > 1) {
    570         for (int i = 0; i < weights->n && nIter > 1; i++) {
    571             psImage *weight = weights->data[i]; // Weight image
    572             if (!weight) {
    573                 // We don't have weights, so no realistic errors: turn off iteration
     568    // 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
    574574                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");
    576576                }
    577577                nIter = 1;
     
    594594                    mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal);
    595595                }
    596                 psImage *weight;        // Weight image
    597                 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];
    599599                } else {
    600600                    errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i];
     
    645645    psFree(images);
    646646    psFree(masks);
    647     psFree(weights);
     647    psFree(variances);
    648648    psFree(refs);
    649649    psFree(regions);
     
    878878        PS_ASSERT_IMAGE_SIZE(readout->mask, data->numCols, data->numRows, NULL);
    879879    }
    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);
    884884    }
    885885
     
    11331133                    mask->data.PS_TYPE_VECTOR_MASK_DATA[r] = (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal);
    11341134                }
    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;
    11371137                } else {
    11381138                    // 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  
    3636#define PEAK_FLUX 1.0e4                 // Peak flux for each source
    3737#define SKY_VALUE 0.0e0                 // Sky value for fake image
    38 #define WEIGHT_VAL 3.0                  // Weighting for image
    39 #define WEIGHT_FACTOR 10.0              // Factor to multiply image by to get weighting
     38#define VARIANCE_VAL 3.0                // Variance for image
     39#define VARIANCE_FACTOR 10.0            // Factor to multiply image by to get variance
    4040#define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF
    4141#define SOURCE_FIT_ITERATIONS 100       // Number of iterations for source fitting
     
    149149            pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source
    150150            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 source
     151            float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source
    152152            if (srcRadius > maxRadius) {
    153153                maxRadius = srcRadius;
     
    216216    psFree(envelope);
    217217
    218     // XXX Setting the weight seems to be an art
     218    // XXX Setting the variance seems to be an art
    219219    // Can't set it too high so that pixels are rejected as insignificant
    220220    // Can't set it too low so that it's hard to get to the minimum
    221221    // 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));
    225226    readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    226227    psImageInit(readout->mask, 0);
     
    242243
    243244        psFree(source->pixels);
    244         psFree(source->weight);
     245        psFree(source->variance);
    245246        psFree(source->maskView);
    246247        psFree(source->maskObj);
    247248        source->pixels = NULL;
    248         source->weight = NULL;
     249        source->variance = NULL;
    249250        source->maskView = NULL;
    250251        source->maskObj = NULL;
     
    280281    options->psfFieldYo = 0;
    281282
    282     pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, WEIGHT_VAL, true);
     283    pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);
    283284
    284285    pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff);
  • branches/pap_branch_20090128/psModules/src/imcombine/pmReadoutCombine.c

    r21183 r21211  
    3737    params->iter = 1;
    3838    params->rej = INFINITY;
    39     params->weights = false;
     39    params->variances = false;
    4040
    4141    return params;
     
    7575    psStringAppend(&comment, "Combining using statistic: %x", params->combine);
    7676    if (!hdu->header) {
    77         hdu->header = psMetadataAlloc();
     77        hdu->header = psMetadataAlloc();
    7878    }
    7979    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     
    8282    // note the clipping parameters, if used
    8383    if (params->combine == PS_STAT_CLIPPED_MEAN) {
    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 weights
    91     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) {
    9292        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", "");
    9494    }
    9595
     
    121121
    122122    // 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);
    124124    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
    125125
     
    177177            return false;
    178178        }
    179         if (params->weights && !readout->weight) {
     179        if (params->variances && !readout->variance) {
    180180            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);
    182182            return false;
    183183        }
     
    190190    }
    191191
    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
    198200    switch (params->combine) {
    199201      case PS_STAT_SAMPLE_MEAN:
     
    250252    psVectorMaskType *maskData = mask->data.PS_TYPE_VECTOR_MASK_DATA;     // Dereference mask
    251253
    252     psVector *weights = NULL;           // Stack of weights
    253     psVector *errors = NULL;            // Stack of errors (sqrt of variance/weights), for psVectorStats
    254     psF32 *weightsData = NULL;          // Dereference weights
    255     if (params->weights) {
    256         weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights
    257         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;
    258260    }
    259261    psVector *index = NULL;             // The indices to sort the pixels
     
    279281    psF32 **outputImage  = output->image->data.F32; // Output image
    280282    psImageMaskType **outputMask   = output->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Output mask
    281     psF32 **outputWeight = NULL; // Output weight map
    282     if (output->weight) {
    283         outputWeight = output->weight->data.F32;
     283    psF32 **outputVariance = NULL; // Output variance map
     284    if (output->variance) {
     285        outputVariance = output->variance->data.F32;
    284286    }
    285287
     
    323325                }
    324326
    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];
    327329                }
    328330
     
    332334                if (scale) {
    333335                    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];
    336338                    }
    337339                }
     
    376378
    377379            // 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");
    380382            }
    381383
     
    386388                outputImage[yOut][xOut] = NAN;
    387389                outputMask[yOut][xOut] = params->blank;
    388                 if (params->weights) {
    389                     outputWeight[yOut][xOut] = NAN;
     390                if (params->variances) {
     391                    outputVariance[yOut][xOut] = NAN;
    390392                }
    391393                sigma->data.F32[yOut][xOut] = NAN;
     
    393395                outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
    394396                outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank;
    395                 if (params->weights) {
     397                if (params->variances) {
    396398                    float stdev = psStatsGetValue(stats, combineStdev);
    397                     outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance
     399                    outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance
    398400                    // XXXX this is not the correct formal error.
    399401                    // also, the weighted mean is not obviously the correct thing here
     
    412414    psFree(pixels);
    413415    psFree(mask);
    414     psFree(weights);
     416    psFree(variances);
    415417    psFree(errors);
    416418    psFree(stats);
  • branches/pap_branch_20090128/psModules/src/imcombine/pmReadoutCombine.h

    r21183 r21211  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2009-01-27 06:39:38 $
     7 * @version $Revision: 1.15.2.1 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-01-29 00:33:51 $
    99 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
    1010 */
     
    2222typedef struct {
    2323    psStatsOptions combine;             ///< Statistic to use when performing the combination
    24     psImageMaskType maskVal;            ///< Mask value
    25     psImageMaskType blank;            ///< Mask value to give blank (i.e., no data) pixels
     24    psImageMaskType maskVal;            ///< Mask value
     25    psImageMaskType blank;            ///< Mask value to give blank (i.e., no data) pixels
    2626    int nKeep;                          ///< Mimimum number of pixels to keep
    2727    float fracHigh;                     ///< Fraction of high pixels to immediately throw
     
    2929    int iter;                           ///< Number of iterations for clipping (for CLIPPED_MEAN only)
    3030    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)?
    3232} pmCombineParams;
    3333
     
    4141/// Combine multiple readouts, applying zero and scale, with optional minmax clipping
    4242bool 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
    4444                      const psVector *zero, ///< Zero corrections to subtract from input, or NULL
    4545                      const psVector *scale, ///< Scale corrections to divide into input, or NULL
  • branches/pap_branch_20090128/psModules/src/imcombine/pmStack.c

    r21183 r21211  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2009-01-27 06:39:38 $
     10 *  @version $Revision: 1.46.2.1 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2009-01-29 00:33:51 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
     
    263263
    264264        psImage *image = data->readout->image; // Image of interest
    265         psImage *variance = data->readout->weight; // Variance ("weight") map of interest
     265        psImage *variance = data->readout->variance; // Variance map of interest
    266266        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
    267267        if (variance) {
     
    442442    *numCols = data->readout->image->numCols;
    443443    *numRows = data->readout->image->numRows;
    444     if (data->readout->weight) {
     444    if (data->readout->variance) {
    445445        *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);
    449449    }
    450450    *haveRejects = (data->reject != NULL);
     
    473473        PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false);
    474474        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);
    478478        }
    479479    }
     
    649649        psImage *combinedImage = combined->image; // Combined image
    650650        psImage *combinedMask = combined->mask; // Combined mask
    651         psImage *combinedVariance = combined->weight; // Combined variance map
     651        psImage *combinedVariance = combined->variance; // Combined variance map
    652652
    653653        psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
     
    702702        }
    703703
    704         psImage *combinedVariance = combined->weight; // Combined variance map
     704        psImage *combinedVariance = combined->variance; // Combined variance map
    705705        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;
    708708        }
    709709
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c

    r21183 r21211  
    6060
    6161    // 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 weight map
     62    // This is required to keep the relative scaling between the image and the variance map
    6363    psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
    6464
     
    287287
    288288// Convolve an image using FFT
    289 static void convolveWeightFFT(psImage *target,// Place the result in here
    290                               psImage *weight, // Weight map to convolve
     289static void convolveVarianceFFT(psImage *target,// Place the result in here
     290                              psImage *variance, // Variance map to convolve
    291291                              psImage *sys, // Systematic error image
    292292                              psImage *mask, // Mask image
     
    302302    bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    303303
    304     psImage *subWeight = convolveSubsetAlloc(weight, border, threaded); // Weight map
     304    psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map
    305305    psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
    306306    psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
    307307
    308308    // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
    309     psImage *convWeight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight
     309    psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
    310310    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
    311311
    312     convolveSubsetFree(weight, subWeight, threaded);
     312    convolveSubsetFree(variance, subVariance, threaded);
    313313    convolveSubsetFree(sys, subSys, threaded);
    314314    convolveSubsetFree(mask, subMask, threaded);
     
    319319        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
    320320            for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
    321                 target->data.F32[yTarget][xTarget] = convWeight->data.F32[ySource][xSource] +
     321                target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] +
    322322                    convSys->data.F32[ySource][xSource];
    323323            }
     
    326326        int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    327327        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
    328             memcpy(&target->data.F32[yTarget][xMin], &convWeight->data.F32[ySource][size], numBytes);
    329         }
    330     }
    331 
    332     psFree(convWeight);
     328            memcpy(&target->data.F32[yTarget][xMin], &convVariance->data.F32[ySource][size], numBytes);
     329        }
     330    }
     331
     332    psFree(convVariance);
    333333    psFree(convSys);
    334334
     
    361361// Convolve a region of an image
    362362static inline void convolveRegion(psImage *convImage, // Convolved image (output)
    363                                   psImage *convWeight, // Convolved weight map (output), or NULL
     363                                  psImage *convVariance, // Convolved variance map (output), or NULL
    364364                                  psImage *convMask, // Convolve mask (output), or NULL
    365365                                  psKernel **kernelImage, // Convolution kernel for the image
    366                                   psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL
     366                                  psKernel **kernelVariance, // Convolution kernel for the variance map, or NULL
    367367                                  psImage *image, // Image to convolve
    368                                   psImage *weight, // Weight map to convolve, or NULL
     368                                  psImage *variance, // Variance map to convolve, or NULL
    369369                                  psImage *sys, // Systematic error image, or NULL
    370370                                  psImage *subMask, // Subtraction mask
     
    381381{
    382382    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
    383     if (weight || subMask) {
    384         *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
     383    if (variance || subMask) {
     384        *kernelVariance = varianceKernel(*kernelVariance, *kernelImage);
    385385    }
    386386
     
    398398    }
    399399
    400     // Convolve the image and weight
     400    // Convolve the image and variance
    401401    if (useFFT) {
    402402        // Use Fast Fourier Transform to do the convolution
    403403        // This provides a big speed-up for large kernels
    404404        convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size);
    405         if (weight) {
    406             convolveWeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size);
     405        if (variance) {
     406            convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size);
    407407        }
    408408    } else {
    409409        // XXX Direct convolution doesn't account for bad pixels yet
    410410        convolveDirect(convImage, image, *kernelImage, region, background, kernels->size);
    411         if (weight) {
    412             convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size);
     411        if (variance) {
     412            convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size);
    413413        }
    414414    }
     
    885885                psFree(stamp->image1);
    886886                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;
    889889                psFree(stamp->matrix1);
    890890                psFree(stamp->matrix2);
     
    989989
    990990
    991 // XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
     991// XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data
    992992static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
    993993                                     int x0, int y0, // Offsets for image
     
    10101010
    10111011    psKernel *kernelImage = NULL;       // Kernel for the images
    1012     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
     1012    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
    10131013
    10141014    // Only generate polynomial values every kernel footprint, since we have already assumed
     
    10201020
    10211021    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,
    10241024                       maskBad, maskPoor, poorFrac, useFFT, false);
    10251025    }
    10261026    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,
    10291029                       maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
    10301030    }
    10311031
    10321032    psFree(kernelImage);
    1033     psFree(kernelWeight);
     1033    psFree(kernelVariance);
    10341034    psFree(polyValues);
    10351035
     
    11481148            }
    11491149        }
    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);
    11531153                if (threaded) {
    1154                     psMutexInit(out1->weight);
     1154                    psMutexInit(out1->variance);
    11551155                }
    11561156            }
    1157             psImageInit(out1->weight, 0.0);
     1157            psImageInit(out1->variance, 0.0);
    11581158        }
    11591159    }
     
    11651165            }
    11661166        }
    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);
    11701170                if (threaded) {
    1171                     psMutexInit(out2->weight);
     1171                    psMutexInit(out2->variance);
    11721172                }
    11731173            }
    1174             psImageInit(out2->weight, 0.0);
     1174            psImageInit(out2->variance, 0.0);
    11751175        }
    11761176    }
     
    12321232    psImage *polyValues = NULL;         // Pre-calculated polynomial values
    12331233    psKernel *kernelImage = NULL;       // Kernel for the images
    1234     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
     1234    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
    12351235#endif
    12361236
     
    13451345        if (out2) {
    13461346            out2->image = psMemIncrRefCounter(ro2->image);
    1347             out2->weight = psMemIncrRefCounter(ro2->weight);
     1347            out2->variance = psMemIncrRefCounter(ro2->variance);
    13481348            out2->mask = psMemIncrRefCounter(ro2->mask);
    13491349        }
     
    13521352        if (out1) {
    13531353            out1->image = psMemIncrRefCounter(ro1->image);
    1354             out1->weight = psMemIncrRefCounter(ro1->weight);
     1354            out1->variance = psMemIncrRefCounter(ro1->variance);
    13551355            out1->mask = psMemIncrRefCounter(ro1->mask);
    13561356        }
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionEquation.c

    r20561 r21211  
    2424static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication
    2525                                         const psKernel *image2, // Second image in multiplication
    26                                          const psKernel *weight, // Weight image
     26                                         const psKernel *variance, // Variance image
    2727                                         int footprint // (Half-)Size of stamp
    2828    )
     
    3131    for (int y = - footprint; y <= footprint; y++) {
    3232        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];
    3434        }
    3535    }
     
    4242                                           const psKernel *image1, // First image in multiplication
    4343                                           const psKernel *image2, // Second image in multiplication
    44                                            const psKernel *weight, // Weight image
     44                                           const psKernel *variance, // Variance image
    4545                                           const psImage *polyValues, // Spatial polynomial values
    4646                                           int numKernels, // Number of kernel basis functions
     
    5050    )
    5151{
    52     double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products
     52    double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products
    5353    if (!isfinite(sum)) {
    5454        return false;
     
    7979                                           const psKernel *image1, // First image in multiplication
    8080                                           const psKernel *image2, // Second image in multiplication
    81                                            const psKernel *weight, // Weight image
     81                                           const psKernel *variance, // Variance image
    8282                                           const psImage *polyValues, // Spatial polynomial values
    8383                                           int numKernels, // Number of kernel basis functions
     
    8787    )
    8888{
    89     double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products
     89    double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products
    9090    if (!isfinite(sum)) {
    9191        return false;
     
    120120                                  const psArray *convolutions1, // Convolutions for element 1
    121121                                  const psArray *convolutions2, // Convolutions for element 2
    122                                   const psKernel *weight, // Weight image
     122                                  const psKernel *variance, // Variance image
    123123                                  const psImage *polyValues, // Polynomial values
    124124                                  int numKernels, // Number of kernel basis functions
     
    135135            psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element
    136136
    137             if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels,
     137            if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels,
    138138                                         footprint, spatialOrder, symmetric)) {
    139139                psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j);
     
    151151                            const psArray *convolutions, // Convolutions of source with kernels
    152152                            const psKernel *input, // Input stamp, or NULL
    153                             const psKernel *weight, // Weight stamp
     153                            const psKernel *variance, // Variance stamp
    154154                            const psImage *polyValues, // Spatial polynomial values
    155155                            int footprint, // (Half-)Size of stamp
     
    171171
    172172    // 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,
    174174                               spatialOrder, footprint)) {
    175175        return false;
     
    185185
    186186            // 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,
    188188                                         footprint, spatialOrder, true)) {
    189189                psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i);
     
    195195            for (int y = - footprint; y <= footprint; y++) {
    196196                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];
    198198                }
    199199            }
     
    218218        for (int y = - footprint; y <= footprint; y++) {
    219219            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];
    221221                double value = input->kernel[y][x] * invNoise2;
    222222                sumI += value;
     
    253253                            const psKernel *input, // Input stamp, or NULL if !normAndBG
    254254                            const psKernel *target, // Target stamp
    255                             const psKernel *weight, // Weight stamp
     255                            const psKernel *variance, // Variance stamp
    256256                            const psImage *polyValues, // Spatial polynomial values
    257257                            int footprint, // (Half-)Size of stamp
     
    277277        for (int y = - footprint; y <= footprint; y++) {
    278278            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];
    280280            }
    281281        }
     
    297297        for (int y = - footprint; y <= footprint; y++) {
    298298            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];
    300300                sumIT += value * input->kernel[y][x];
    301301                sumT += value;
     
    329329                                 const psArray *convolutions2, // Convolutions of image 2
    330330                                 const psKernel *image1, // Image 1 stamp
    331                                  const psKernel *weight, // Weight stamp
     331                                 const psKernel *variance, // Variance stamp
    332332                                 const psImage *polyValues, // Spatial polynomial values
    333333                                 int footprint // (Half-)Size of stamp
     
    348348    PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels);
    349349
    350     if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels,
     350    if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels,
    351351                               spatialOrder, footprint)) {
    352352        return false;
     
    356356        // Normalisation
    357357        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,
    359359                                     footprint, spatialOrder, false)) {
    360360            psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i);
     
    366366        for (int y = - footprint; y <= footprint; y++) {
    367367            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];
    369369            }
    370370        }
     
    559559      case PM_SUBTRACTION_MODE_1:
    560560        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
    561                                  stamp->weight, polyValues, footprint, true);
     561                                 stamp->variance, polyValues, footprint, true);
    562562        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);
    564564        break;
    565565      case PM_SUBTRACTION_MODE_2:
    566566        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
    567                                  stamp->weight, polyValues, footprint, true);
     567                                 stamp->variance, polyValues, footprint, true);
    568568        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);
    570570        break;
    571571      case PM_SUBTRACTION_MODE_DUAL:
     
    581581#endif
    582582        status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
    583                                   stamp->weight, polyValues, footprint, true);
     583                                  stamp->variance, polyValues, footprint, true);
    584584        status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
    585                                   stamp->weight, polyValues, footprint, false);
     585                                  stamp->variance, polyValues, footprint, false);
    586586        status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
    587                                        stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
     587                                       stamp->convolutions2, stamp->image1, stamp->variance, polyValues,
    588588                                       footprint);
    589589        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);
    591591        status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
    592                                   stamp->image2, stamp->weight, polyValues, footprint, false);
     592                                  stamp->image2, stamp->variance, polyValues, footprint, false);
    593593        break;
    594594      default:
     
    10331033
    10341034        // Calculate residuals
    1035         psKernel *weight = stamp->weight; // Weight postage stamp
     1035        psKernel *variance = stamp->variance; // Variance postage stamp
    10361036        psImageInit(residual->image, 0.0);
    10371037        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     
    10981098        for (int y = - footprint; y <= footprint; y++) {
    10991099            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];
    11011101                deviation += dev;
    11021102#ifdef TESTING
     
    11411141            psFitsClose(fits);
    11421142        }
    1143         if (stamp->weight) {
     1143        if (stamp->variance) {
    11441144            psString filename = NULL;
    1145             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
     1145            psStringAppend(&filename, "stamp_variance_%03d.fits", i);
    11461146            psFits *fits = psFitsOpen(filename, "w");
    11471147            psFree(filename);
    1148             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
     1148            psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL);
    11491149            psFitsClose(fits);
    11501150        }
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionMatch.c

    r21183 r21211  
    6060                      const pmReadout *ro2, // Readout 2
    6161                      const psImage *subMask, // Mask for subtraction, or NULL
    62                       psImage *weight,  // Weight map
     62                      psImage *variance,  // Variance map
    6363                      const psRegion *region, // Region of interest, or NULL
    6464                      float threshold,  // Threshold for stamp finding
     
    8080
    8181    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)) {
    8383        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    8484        return false;
     
    111111            conv1->mask = NULL;
    112112        }
    113         if (conv1->weight) {
    114             psFree(conv1->weight);
    115             conv1->weight = NULL;
     113        if (conv1->variance) {
     114            psFree(conv1->variance);
     115            conv1->variance = NULL;
    116116        }
    117117    }
     
    126126            conv2->mask = NULL;
    127127        }
    128         if (conv2->weight) {
    129             psFree(conv2->weight);
    130             conv2->weight = NULL;
     128        if (conv2->variance) {
     129            psFree(conv2->variance);
     130            conv2->variance = NULL;
    131131        }
    132132    }
     
    191191    }
    192192
    193     // Where does our weight map come from?
    194     // Getting the weight exactly right is not necessary --- it's just used for weighting.
    195     psImage *weight = NULL;             // Weight image to use
    196     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);
    202202    } else {
    203         weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
     203        variance = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
    204204    }
    205205
     
    274274            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    275275            // 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,
    277277                           size, footprint, subMode)) {
    278278                goto MATCH_ERROR;
     
    336336                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    337337
    338                 if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
     338                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing,
    339339                               size, footprint, subMode)) {
    340340                    goto MATCH_ERROR;
     
    439439    psFree(subMask);
    440440    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)) {
    445445        psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
    446446        goto MATCH_ERROR;
     
    482482    psFree(kernels);
    483483    psFree(stamps);
    484     psFree(weight);
     484    psFree(variance);
    485485    psFree(rng);
    486486    return false;
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionParams.c

    r18287 r21211  
    7171                            double *sumII, // Sum of I(x)^2/sigma(x)^2
    7272                            double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
    73                             const pmSubtractionStamp *stamp, // Stamp with weight
     73                            const pmSubtractionStamp *stamp, // Stamp with variance
    7474                            const psKernel *target, // Target stamp
    7575                            int kernelIndex, // Index for kernel component
     
    7878    )
    7979{
    80     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
     80    psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
    8181    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    8282
    8383    for (int y = -footprint; y <= footprint; y++) {
    8484        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    85         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     85        psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
    8686        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    8787        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
     
    9898static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2
    9999                                   double *sumCC, // Sum of conv(x)^2/sigma(x)^2
    100                                    const pmSubtractionStamp *stamp, // Stamp with input and weight
     100                                   const pmSubtractionStamp *stamp, // Stamp with input and variance
    101101                                   int kernelIndex, // Index for kernel component
    102102                                   int footprint, // Size of region of interest
     
    104104    )
    105105{
    106     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
     106    psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
    107107    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    108108
    109109    for (int y = -footprint; y <= footprint; y++) {
    110         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     110        psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
    111111        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    112112        for (int x = -footprint; x <= footprint; x++, wt++, conv++) {
     
    120120
    121121static double accumulateChi2(const psKernel *target, // Target stamp
    122                              pmSubtractionStamp *stamp, // Stamp with weight
     122                             pmSubtractionStamp *stamp, // Stamp with variance
    123123                             int kernelIndex, // Index for kernel component
    124124                             double coeff, // Coefficient of convolution
     
    129129{
    130130    double chi2 = 0.0;
    131     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
     131    psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
    132132    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    133133
    134134    for (int y = -footprint; y <= footprint; y++) {
    135135        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    136         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     136        psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
    137137        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    138138        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
     
    146146// Return the initial value of chi^2
    147147static double initialChi2(const psKernel *target, // Target stamp
    148                           const pmSubtractionStamp *stamp, // Stamp with weight
     148                          const pmSubtractionStamp *stamp, // Stamp with variance
    149149                          int footprint, // Size of convolution
    150150                          pmSubtractionMode mode // Mode of subtraction
    151151    )
    152152{
    153     psKernel *weight = stamp->weight;   // Weight map
     153    psKernel *variance = stamp->variance;   // Variance map
    154154    psKernel *source;                   // Source stamp
    155155    switch (mode) {
     
    167167    for (int y = -footprint; y <= footprint; y++) {
    168168        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    169         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     169        psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
    170170        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
    171171        for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
     
    180180// Subtract a convolution from the input
    181181static void subtractConvolution(psKernel *target, // Target stamp
    182                                 const pmSubtractionStamp *stamp, // Stamp with weight
     182                                const pmSubtractionStamp *stamp, // Stamp with variance
    183183                                int kernelIndex, // Index for kernel component
    184184                                float coeff, // Coefficient of subtraction
     
    288288
    289289        // This sum is invariant to the kernel
    290         psKernel *weight = stamp->weight; // Weight map for stamp
     290        psKernel *variance = stamp->variance; // Variance map for stamp
    291291        for (int v = -footprint; v <= footprint; v++) {
    292             psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map
     292            psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map
    293293            for (int u = -footprint; u <= footprint; u++, wt++) {
    294294                sum1 += 1.0 / *wt;
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionStamps.c

    r21183 r21211  
    5454    psFree(stamp->image1);
    5555    psFree(stamp->image2);
    56     psFree(stamp->weight);
     56    psFree(stamp->variance);
    5757    psFree(stamp->convolutions1);
    5858    psFree(stamp->convolutions2);
     
    211211    stamp->image1 = NULL;
    212212    stamp->image2 = NULL;
    213     stamp->weight = NULL;
     213    stamp->variance = NULL;
    214214    stamp->convolutions1 = NULL;
    215215    stamp->convolutions2 = NULL;
     
    333333                psFree(stamp->image1);
    334334                psFree(stamp->image2);
    335                 psFree(stamp->weight);
     335                psFree(stamp->variance);
    336336                psFree(stamp->convolutions1);
    337337                psFree(stamp->convolutions2);
    338                 stamp->image1 = stamp->image2 = stamp->weight = NULL;
     338                stamp->image1 = stamp->image2 = stamp->variance = NULL;
    339339                stamp->convolutions1 = stamp->convolutions2 = NULL;
    340340
     
    480480
    481481bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
    482                                 psImage *weight, int kernelSize)
     482                                psImage *variance, int kernelSize)
    483483{
    484484    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     
    490490        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
    491491    }
    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);
    495495    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    496496
     
    520520        assert(stamp->image1 == NULL);
    521521        assert(stamp->image2 == NULL);
    522         assert(stamp->weight == NULL);
     522        assert(stamp->variance == NULL);
    523523
    524524        psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
     
    534534        }
    535535
    536         psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
    537         stamp->weight = psKernelAllocFromImage(wtSub, size, size);
     536        psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp
     537        stamp->variance = psKernelAllocFromImage(wtSub, size, size);
    538538        psFree(wtSub);                  // Drop reference
    539539
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionStamps.h

    r20465 r21211  
    6060    psKernel *image1;                   ///< Reference image postage stamp
    6161    psKernel *image2;                   ///< Input image postage stamp
    62     psKernel *weight;                   ///< Weight image postage stamp, or NULL
     62    psKernel *variance;                 ///< Variance image postage stamp, or NULL
    6363    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
    6464    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
     
    120120                                psImage *image1, ///< Reference image
    121121                                psImage *image2, ///< Input image (or NULL)
    122                                 psImage *weight, ///< Weight (variance) map
     122                                psImage *variance, ///< Variance map
    123123                                int kernelSize ///< Kernel half-size
    124124    );
  • branches/pap_branch_20090128/psModules/src/imcombine/pmSubtractionThreads.c

    r19765 r21211  
    2727        psMutexInit(ro->image);
    2828    }
    29     if (ro->weight) {
    30         psMutexInit(ro->weight);
     29    if (ro->variance) {
     30        psMutexInit(ro->variance);
    3131    }
    3232
     
    4343        psMutexDestroy(ro->image);
    4444    }
    45     if (ro->weight) {
    46         psMutexDestroy(ro->weight);
     45    if (ro->variance) {
     46        psMutexDestroy(ro->variance);
    4747    }
    4848
  • branches/pap_branch_20090128/psModules/src/objects/pmSource.c

    r21183 r21211  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.67.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4444    psFree(tmp->peak);
    4545    psFree(tmp->pixels);
    46     psFree(tmp->weight);
     46    psFree(tmp->variance);
    4747    psFree(tmp->maskObj);
    4848    psFree(tmp->maskView);
     
    6666
    6767    psFree (source->pixels);
    68     psFree (source->weight);
     68    psFree (source->variance);
    6969    psFree (source->maskObj);
    7070    psFree (source->maskView);
     
    7373
    7474    source->pixels = NULL;
    75     source->weight = NULL;
     75    source->variance = NULL;
    7676    source->maskObj = NULL;
    7777    source->maskView = NULL;
     
    101101    source->peak = NULL;
    102102    source->pixels = NULL;
    103     source->weight = NULL;
     103    source->variance = NULL;
    104104    source->maskObj = NULL;
    105105    source->maskView = NULL;
     
    144144    }
    145145    // this copy is used to allow multiple fit attempts on the
    146     // same object.  the pixels, weight, and mask arrays all point to
     146    // same object.  the pixels, variance, and mask arrays all point to
    147147    // the same original subarrays.  the peak and moments point at
    148148    // the original values.
     
    167167    // We want a new view, but pointing at the same pixels.
    168168    source->pixels   = psImageCopyView(NULL, in->pixels);
    169     source->weight   = psImageCopyView(NULL, in->weight);
     169    source->variance   = psImageCopyView(NULL, in->variance);
    170170    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
    171171
     
    199199    // these images are subset images of the equivalent parents
    200200    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);
    203203    }
    204204    if (readout->mask) {
     
    237237
    238238    extend |= (mySource->pixels == NULL);
    239     extend |= (mySource->weight == NULL);
     239    extend |= (mySource->variance == NULL);
    240240    extend |= (mySource->maskObj == NULL);
    241241    extend |= (mySource->maskView == NULL);
     
    245245        // re-create the subimage
    246246        psFree (mySource->pixels);
    247         psFree (mySource->weight);
     247        psFree (mySource->variance);
    248248        psFree (mySource->maskView);
    249249
    250250        mySource->pixels   = psImageSubset(readout->image,  newRegion);
    251         mySource->weight   = psImageSubset(readout->weight, newRegion);
     251        mySource->variance   = psImageSubset(readout->variance, newRegion);
    252252        mySource->maskView = psImageSubset(readout->mask,   newRegion);
    253253        mySource->region   = newRegion;
     
    670670    pmSource->peak
    671671    pmSource->pixels
    672     pmSource->weight
     672    pmSource->variance
    673673    pmSource->mask
    674674
     
    737737
    738738        psF32 *vPix = source->pixels->data.F32[row];
    739         psF32 *vWgt = source->weight->data.F32[row];
     739        psF32 *vWgt = source->variance->data.F32[row];
    740740        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    741741
     
    923923        if (mode & PM_MODEL_OP_NOISE) {
    924924            // XXX need to scale by the gain...
    925             target = source->weight->data.F32;
     925            target = source->variance->data.F32;
    926926        }
    927927
     
    945945    psImage *target = source->pixels;
    946946    if (mode & PM_MODEL_OP_NOISE) {
    947         target = source->weight;
     947        target = source->variance;
    948948    }
    949949
  • branches/pap_branch_20090128/psModules/src/objects/pmSource.h

    r21183 r21211  
    33 * @author EAM, IfA; GLG, MHPCC
    44 *
    5  * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
    6  * @date $Date: 2009-01-27 06:39:38 $
     5 * @version $Revision: 1.27.2.1 $ $Name: not supported by cvs2svn $
     6 * @date $Date: 2009-01-29 00:33:51 $
    77 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    88 */
     
    7272    pmPeak  *peak;                      ///< Description of peak pixel.
    7373    psImage *pixels;                    ///< Rectangular region including object pixels.
    74     psImage *weight;                    ///< Image variance.
     74    psImage *variance;                    ///< Image variance.
    7575    psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
    7676    psImage *maskView;                  ///< view into global image mask for this object region
     
    210210    psMetadata *metadata,               ///< Contains classification parameters
    211211    pmPSFClump clump,                   ///< Statistics about the PSF clump
    212     psImageMaskType maskSat             ///< Mask value for saturated pixels
     212    psImageMaskType maskSat             ///< Mask value for saturated pixels
    213213);
    214214
  • branches/pap_branch_20090128/psModules/src/objects/pmSourceFitModel.c

    r21183 r21211  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.29.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    6262    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    6363    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);
    6565
    6666    psBool fitStatus = true;
     
    8484                continue;
    8585            }
    86             // skip zero-weight points
    87             if (source->weight->data.F32[i][j] == 0) {
     86            // skip zero-variance points
     87            if (source->variance->data.F32[i][j] == 0) {
    8888                continue;
    8989            }
     
    102102
    103103            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
    104             // as weight to avoid the bias from systematic errors here we would just use the
     104            // as variance to avoid the bias from systematic errors here we would just use the
    105105            // source sky variance
    106106            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];
    108108            } else {
    109109                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
  • branches/pap_branch_20090128/psModules/src/objects/pmSourceFitSet.c

    r21183 r21211  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.14.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 *
     
    5555
    5656    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
    6262    // before pmSourceFitSetInit is called
    6363    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?");
    6565    }
    6666    return true;
     
    125125    // we do not need to lock to do this....
    126126    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;
    133133    }
    134134    psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared");
     
    137137
    138138    pthread_mutex_lock (&fitSetInitMutex);
    139          
     139
    140140    // find an entry that is NULL and place it there
    141141    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;
    146146    }
    147147    pthread_mutex_unlock (&fitSetInitMutex);
     
    160160    pmSourceFitSetData *thisSet = NULL;
    161161    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;
    168168    }
    169169    psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found");
     
    184184    pmSourceFitSetData *thisSet = NULL;
    185185    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;
    192192    }
    193193    psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found");
     
    256256        psVector *derivOne = set->derivSet->data[i];
    257257
    258         // 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
     258        // 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
    262262        assert ((deriv == NULL) || (paramOne->n == derivOne->n));
    263        
     263
    264264        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            }
    268268            if (deriv) {
    269269                deriv->data.F32[n] = derivOne->data.F32[j];
     
    275275
    276276// distribute parameters from single param and deriv vectors into FitSet models
    277 bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) 
     277bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param)
    278278{
    279279    PS_ASSERT_VECTOR_NON_NULL(param, false);
     
    301301
    302302// set the model parameters for this fit set
    303 bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, 
     303bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam,
    304304                           const psVector *param, pmSource *source,
    305305                           psMinimization *myMin, int nPix, bool fitStatus)
     
    420420
    421421    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        }
    427427    }
    428428    return true;
     
    435435        psVector *derivOne = set->derivSet->data[i];
    436436        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]);
    438438        }
    439439    }
     
    450450    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    451451    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);
    453453
    454454    bool fitStatus = true;
     
    471471                continue;
    472472            }
    473             // skip zero-weight points
    474             if (source->weight->data.F32[i][j] == 0) {
     473            // skip zero-variance points
     474            if (source->variance->data.F32[i][j] == 0) {
    475475                continue;
    476476            }
     
    489489
    490490            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
    491             // as weight to avoid the bias from systematic errors here we would just use the
     491            // as variance to avoid the bias from systematic errors here we would just use the
    492492            // source sky variance
    493493            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];
    495495            } else {
    496496                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
     
    543543        psFree (params);
    544544        psFree(constraint);
    545         pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
     545        pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
    546546        return(false);
    547547    }
  • branches/pap_branch_20090128/psModules/src/objects/pmSourceMoments.c

    r21183 r21211  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.7.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4747    pmSource->peak
    4848    pmSource->pixels
    49     pmSource->weight (optional)
     49    pmSource->variance (optional)
    5050    pmSource->mask (optional)
    5151
     
    9999
    100100        psF32 *vPix = source->pixels->data.F32[row];
    101         psF32 *vWgt = source->weight->data.F32[row];
     101        psF32 *vWgt = source->variance->data.F32[row];
    102102        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    103103
     
    110110                vMsk++;
    111111            }
    112             if (isnan(*vPix)) continue;
     112            if (isnan(*vPix)) continue;
    113113
    114114            psF32 xDiff = col - xPeak;
     
    189189
    190190        psF32 *vPix = source->pixels->data.F32[row];
    191         psF32 *vWgt = source->weight->data.F32[row];
     191        psF32 *vWgt = source->variance->data.F32[row];
    192192        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    193193
     
    200200                vMsk++;
    201201            }
    202             if (isnan(*vPix)) continue;
     202            if (isnan(*vPix)) continue;
    203203
    204204            psF32 xDiff = col - xCM;
     
    206206
    207207            // 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);
    210210            if (r > radius) continue;
    211211
     
    278278
    279279    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,
    281281             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    282282             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     
    286286
    287287    if (source->moments->Mxx < 0) {
    288         fprintf (stderr, "error: neg second moment??\n");
     288        fprintf (stderr, "error: neg second moment??\n");
    289289    }
    290290    if (source->moments->Myy < 0) {
    291         fprintf (stderr, "error: neg second moment??\n");
     291        fprintf (stderr, "error: neg second moment??\n");
    292292    }
    293293
     
    341341
    342342        psF32 *vPix = source->pixels->data.F32[row];
    343         psF32 *vWgt = source->weight->data.F32[row];
     343        psF32 *vWgt = source->variance->data.F32[row];
    344344        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    345345
     
    352352                vMsk++;
    353353            }
    354             if (isnan(*vPix)) continue;
     354            if (isnan(*vPix)) continue;
    355355
    356356            psF32 xDiff = col - xPeak;
     
    407407    }
    408408
    409 # if (PS_TRACE_ON) 
     409# if (PS_TRACE_ON)
    410410    float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0);
    411411    float Sxy = XY / Sum;
  • branches/pap_branch_20090128/psModules/src/objects/pmSourcePhotometry.c

    r21183 r21211  
    33 *  @author EAM, IfA; GLG, MHPCC
    44 *
    5  *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2009-01-27 06:39:38 $
     5 *  @version $Revision: 1.48.2.1 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2009-01-29 00:33:51 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    446446    assert (Pj != NULL);
    447447
    448     const psImage *Wi = Mi->weight;
     448    const psImage *Wi = Mi->variance;
    449449    if (!unweighted_sum) {
    450450        assert (Wi != NULL);
     
    511511    assert (Pj != NULL);
    512512
    513     const psImage *Wi = Mi->weight;
     513    const psImage *Wi = Mi->variance;
    514514    if (!unweighted_sum) {
    515515        assert (Wi != NULL);
     
    567567    const psImage *Pi = Mi->pixels;
    568568    assert (Pi != NULL);
    569     const psImage *Wi = Mi->weight;
     569    const psImage *Wi = Mi->variance;
    570570    if (!unweighted_sum) {
    571571        assert (Wi != NULL);
     
    612612# endif
    613613
    614 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight,
     614bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance,
    615615                    psImageMaskType maskVal)
    616616{
     
    618618    PS_ASSERT_PTR_NON_NULL(image, false);
    619619    PS_ASSERT_PTR_NON_NULL(mask, false);
    620     PS_ASSERT_PTR_NON_NULL(weight, false);
     620    PS_ASSERT_PTR_NON_NULL(variance, false);
    621621
    622622    double dC = 0.0;
     
    626626            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)
    627627                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];
    631631            Npix ++;
    632632        }
     
    648648    const psImage *Pi = Mi->modelFlux;
    649649    assert (Pi != NULL);
    650     const psImage *Wi = Mi->weight;
     650    const psImage *Wi = Mi->variance;
    651651    if (!unweighted_sum) {
    652652        assert (Wi != NULL);
     
    706706    assert (Pj != NULL);
    707707
    708     const psImage *Wi = Mi->weight;
     708    const psImage *Wi = Mi->variance;
    709709    if (!unweighted_sum) {
    710710        assert (Wi != NULL);
     
    770770    assert (Pj != NULL);
    771771
    772     const psImage *Wi = Mi->weight;
     772    const psImage *Wi = Mi->variance;
    773773    if (!unweighted_sum) {
    774774        assert (Wi != NULL);
  • branches/pap_branch_20090128/psModules/src/objects/pmSourceSky.c

    r21183 r21211  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-01-27 06:39:38 $
     8 *  @version $Revision: 1.19.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-01-29 00:33:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    113113    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
    114114    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);
    116116    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
    117117    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     
    127127    }
    128128
    129     psImage *image = source->weight;
     129    psImage *image = source->variance;
    130130    psImage *mask  = source->maskObj;
    131131    pmPeak *peak  = source->peak;
Note: See TracChangeset for help on using the changeset viewer.