IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21363


Ignore:
Timestamp:
Feb 5, 2009, 4:31:25 PM (17 years ago)
Author:
Paul Price
Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

Location:
trunk/psModules/src
Files:
47 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/camera/pmFPA.c

    r19013 r21363  
    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;

  • trunk/psModules/src/camera/pmFPA.h

    r19013 r21363  
    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.26 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2009-02-06 02:31:24 $
    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/// @}

  • trunk/psModules/src/camera/pmFPAConstruct.c

    r17911 r21363  
    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

  • trunk/psModules/src/camera/pmFPACopy.c

    r20634 r21363  
    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)) {

  • trunk/psModules/src/camera/pmFPAExtent.c

    r20635 r21363  
    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

  • trunk/psModules/src/camera/pmFPAMaskWeight.c

    r21183 r21363  
    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

  • trunk/psModules/src/camera/pmFPAMaskWeight.h

    r21183 r21363  
    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.18 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-02-06 02:31:24 $
    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

  • trunk/psModules/src/camera/pmFPAMosaic.c

    r21183 r21363  
    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

  • trunk/psModules/src/camera/pmFPARead.c

    r21183 r21363  
    99#include <pslib.h>
    1010
     11#include "pmConfig.h"
     12#include "pmConfigMask.h"
    1113#include "pmHDU.h"
    1214#include "pmFPA.h"
    … …  
    2830    FPA_READ_TYPE_IMAGE,                // Read image
    2931    FPA_READ_TYPE_MASK,                 // Read mask
    30     FPA_READ_TYPE_WEIGHT,               // Read weight map
     32    FPA_READ_TYPE_VARIANCE,             // Read variance map
    3133    FPA_READ_TYPE_HEADER                // Read header
    3234} fpaReadType;
    … …  
    5456      case FPA_READ_TYPE_MASK:
    5557        return readout->thisMaskScan;
    56       case FPA_READ_TYPE_WEIGHT:
    57         return readout->thisWeightScan;
     58      case FPA_READ_TYPE_VARIANCE:
     59        return readout->thisVarianceScan;
    5860      default:
    5961        psAbort("Unknown read type: %x\n", type);
    … …  
    7476        readout->thisMaskScan = thisScan;
    7577        return readout->lastMaskScan;
    76       case FPA_READ_TYPE_WEIGHT:
    77         readout->thisWeightScan = thisScan;
    78         return readout->lastWeightScan;
     78      case FPA_READ_TYPE_VARIANCE:
     79        readout->thisVarianceScan = thisScan;
     80        return readout->lastVarianceScan;
    7981      default:
    8082        psAbort("Unknown read type: %x\n", type);
    … …  
    9395      case FPA_READ_TYPE_MASK:
    9496        return readout->lastMaskScan;
    95       case FPA_READ_TYPE_WEIGHT:
    96         return readout->lastWeightScan;
     97      case FPA_READ_TYPE_VARIANCE:
     98        return readout->lastVarianceScan;
    9799      default:
    98100        psAbort("Unknown read type: %x\n", type);
    … …  
    113115        readout->lastMaskScan = lastScan;
    114116        return readout->lastMaskScan;
    115       case FPA_READ_TYPE_WEIGHT:
    116         readout->lastWeightScan = lastScan;
    117         return readout->lastWeightScan;
     117      case FPA_READ_TYPE_VARIANCE:
     118        readout->lastVarianceScan = lastScan;
     119        return readout->lastVarianceScan;
    118120      default:
    119121        psAbort("Unknown read type: %x\n", type);
    … …  
    132134      case FPA_READ_TYPE_MASK:
    133135        return &readout->mask;
    134       case FPA_READ_TYPE_WEIGHT:
    135         return &readout->weight;
     136      case FPA_READ_TYPE_VARIANCE:
     137        return &readout->variance;
    136138      default:
    137139        psAbort("Unknown read type: %x\n", type);
    … …  
    193195
    194196    return naxis3;
     197}
     198
     199// Determine whether a FITS file contains covariance matrices
     200static bool hduCovariance(pmHDU *hdu,   // Header data unit
     201                          psFits *fits  // FITS file
     202    )
     203{
     204    if (hdu->extname && !psFitsMoveExtName(fits, hdu->extname)) {
     205        psError(PS_ERR_IO, false, "Unable to move to extension %s", hdu->extname);
     206        return false;
     207    }
     208    // Need to explicitly read the header, since the HDU may not contain the variance header
     209    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
     210    if (!header) {
     211        psError(PS_ERR_IO, false, "Unable to read variance header.");
     212        return false;
     213    }
     214    bool mdok;                          // Status of MD lookup
     215    bool covar = psMetadataLookupBool(&mdok, header, PM_HDU_COVARIANCE_KEYWORD); // Got covariance?
     216    psFree(header);
     217    return covar;
    195218}
    196219
    … …  
    350373    *target = psImageSubset(image, region);
    351374
    352     // Get the list of overscans: only for IMAGE types (no overscan for MASK and WEIGHT)
     375    // Get the list of overscans: only for IMAGE types (no overscan for MASK and VARIANCE)
    353376    if (type == FPA_READ_TYPE_IMAGE) {
    354377        if (readout->bias->n != 0) {
    … …  
    515538    }
    516539
    517     // XXX for IMAGE, we need the CELL.BAD value, but for MASK, we need the BAD mask value
    518 
    519     float bad = 0;
    520     if (type == FPA_READ_TYPE_MASK) {
    521       bad = 1.0;
    522     } else {
    523       bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
     540    // Need to set the invalid (unread) pixels appropriately, and to read the mask bits
     541    float bad = 0;                      // Bad level
     542    switch (type) {
     543      case FPA_READ_TYPE_MASK: {
     544          // Need to explicitly read the header, since what's in the pmHDU may not correspond to the mask
     545          psMetadata *header = psFitsReadHeader(NULL, fits);
     546          if (!header) {
     547              psError(PS_ERR_IO, false, "Unable to read mask header.");
     548              return false;
     549          }
     550          if (!pmConfigMaskReadHeader(config, header)) {
     551              psError(PS_ERR_IO, false, "Unable to determine mask bits");
     552              psFree(header);
     553              return false;
     554          }
     555          psFree(header);
     556          bad = pmConfigMaskGet("BAD", config);
     557          break;
     558      }
     559      case FPA_READ_TYPE_IMAGE:
     560      case FPA_READ_TYPE_VARIANCE:
     561        bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD");
     562        break;
     563      default:
     564        psAbort("Unrecognised type: %x", type);
    524565    }
    525566
    … …  
    582623    *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]);
    583624
    584     // Read overscans only for "image" type --- weights and masks shouldn't record overscans
     625    // Read overscans only for "image" type --- variances and masks shouldn't record overscans
    585626    if (type == FPA_READ_TYPE_IMAGE) {
    586627        // Blow away existing data
    … …  
    608649}
    609650
    610 // Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadWeight
     651// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadVariance
    611652// Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts.
    612653static bool cellRead(pmCell *cell,      // Cell into which to read
    … …  
    640681        dataPointer = hdu->masks;
    641682        break;
    642       case FPA_READ_TYPE_WEIGHT:
    643         hduReadFunc = pmHDUReadWeight;
    644         dataPointer = hdu->weights;
     683      case FPA_READ_TYPE_VARIANCE:
     684        hduReadFunc = pmHDUReadVariance;
     685        dataPointer = hdu->variances;
    645686        break;
    646687      default:
    … …  
    680721        imageArray = hdu->masks;
    681722        break;
    682       case FPA_READ_TYPE_WEIGHT:
    683         imageArray = hdu->weights;
     723      case FPA_READ_TYPE_VARIANCE:
     724        imageArray = hdu->variances;
    684725        break;
    685726      default:
    … …  
    729770
    730771
    731 // Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadWeight
     772// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadVariance
    732773// Iterates over component cells, reading each
    733774static bool chipRead(pmChip *chip,      // Chip into which to read
    … …  
    760801
    761802
    762 // Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadWeight
     803// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadVariance
    763804// Iterates over component chips, reading each
    764805static bool fpaRead(pmFPA *fpa,         // FPA into which to read
    … …  
    852893    float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
    853894    if (!mdok) {
    854         psLogMsg(__func__, PS_LOG_WARN, "CELL.BAD is not set --- assuming zero.\n");
     895        psWarning("CELL.BAD is not set --- assuming zero.\n");
    855896        bad = 0.0;
    856897    }
    … …  
    10911132
    10921133//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1093 // Reading the weight map
    1094 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1095 
    1096 bool pmReadoutMoreWeight(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
     1134// Reading the variance map
     1135//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     1136
     1137bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
    10971138{
    10981139    PS_ASSERT_PTR_NON_NULL(readout, false);
    10991140    PS_ASSERT_FITS_NON_NULL(fits, false);
    11001141
    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,
     1142    return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config);
     1143}
     1144
     1145bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
    11051146                              pmConfig *config)
    11061147{
    … …  
    11101151    PS_ASSERT_INT_NONNEGATIVE(numScans, false);
    11111152
    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)
     1153    return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config);
     1154}
     1155
     1156bool pmReadoutReadVariance(pmReadout *readout, psFits *fits, int z, pmConfig *config)
    11161157{
    11171158    PS_ASSERT_PTR_NON_NULL(readout, false);
    11181159    PS_ASSERT_FITS_NON_NULL(fits, false);
    11191160
    1120     return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_WEIGHT, config);
    1121 }
    1122 
    1123 bool pmCellReadWeight(pmCell *cell, psFits *fits, pmConfig *config)
     1161    return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config);
     1162}
     1163
     1164bool pmCellReadVariance(pmCell *cell, psFits *fits, pmConfig *config)
    11241165{
    11251166    PS_ASSERT_PTR_NON_NULL(cell, false);
    11261167    PS_ASSERT_FITS_NON_NULL(fits, false);
    11271168
    1128     return cellRead(cell, fits, config, FPA_READ_TYPE_WEIGHT);
    1129 }
    1130 
    1131 bool pmChipReadWeight(pmChip *chip, psFits *fits, pmConfig *config)
     1169    if (!cellRead(cell, fits, config, FPA_READ_TYPE_VARIANCE)) {
     1170        return false;
     1171    }
     1172    pmHDU *hdu = pmHDUFromCell(cell);   // Header data unit
     1173    if (hduCovariance(hdu, fits)) {
     1174        return pmCellReadCovariance(cell, fits);
     1175    }
     1176    return true;
     1177}
     1178
     1179bool pmChipReadVariance(pmChip *chip, psFits *fits, pmConfig *config)
    11321180{
    11331181    PS_ASSERT_PTR_NON_NULL(chip, false);
    11341182    PS_ASSERT_FITS_NON_NULL(fits, false);
    11351183
    1136     return chipRead(chip, fits, config, FPA_READ_TYPE_WEIGHT);
    1137 }
    1138 
    1139 bool pmFPAReadWeight(pmFPA *fpa, psFits *fits, pmConfig *config)
     1184    if (!chipRead(chip, fits, config, FPA_READ_TYPE_VARIANCE)) {
     1185        return false;
     1186    }
     1187    pmHDU *hdu = pmHDUFromChip(chip);   // Header data unit
     1188    if (hduCovariance(hdu, fits)) {
     1189        return pmChipReadCovariance(chip, fits);
     1190    }
     1191    return true;
     1192}
     1193
     1194bool pmFPAReadVariance(pmFPA *fpa, psFits *fits, pmConfig *config)
    11401195{
    11411196    PS_ASSERT_PTR_NON_NULL(fpa, false);
    11421197    PS_ASSERT_FITS_NON_NULL(fits, false);
    11431198
    1144     return fpaRead(fpa, fits, config, FPA_READ_TYPE_WEIGHT);
     1199    if (!fpaRead(fpa, fits, config, FPA_READ_TYPE_VARIANCE)) {
     1200        return false;
     1201    }
     1202    pmHDU *hdu = pmHDUFromFPA(fpa);     // Header data unit
     1203    if (hduCovariance(hdu, fits)) {
     1204        return pmFPAReadCovariance(fpa, fits);
     1205    }
     1206    return true;
    11451207}
    11461208
    … …  
    12681330    return numRead;
    12691331}
     1332
     1333
     1334//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     1335// Reading covariance matrices
     1336//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     1337
     1338bool pmCellReadCovariance(pmCell *cell, psFits *fits)
     1339{
     1340    PS_ASSERT_PTR_NON_NULL(cell, false);
     1341    PS_ASSERT_FITS_NON_NULL(fits, false);
     1342
     1343    const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip
     1344    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
     1345    psString extname = NULL;            // Extension name
     1346    psStringAppend(&extname, "COVAR_%s_%s", chipName, cellName);
     1347
     1348    if (!psFitsMoveExtName(fits, extname)) {
     1349        psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
     1350        psFree(extname);
     1351        return false;
     1352    }
     1353    psFree(extname);
     1354
     1355    psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header
     1356    if (!header) {
     1357        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname);
     1358        psFree(header);
     1359        return false;
     1360    }
     1361
     1362    bool mdok;                          // Status of MD lookup
     1363    int x0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.X"); // Centre of matrix in x
     1364    if (!mdok) {
     1365        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre");
     1366        psFree(header);
     1367        return false;
     1368    }
     1369    int y0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.Y"); // Centre of matrix in y
     1370    if (!mdok) {
     1371        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre");
     1372        psFree(header);
     1373        return false;
     1374    }
     1375    psFree(header);
     1376
     1377    psArray *images = psFitsReadImageCube(fits, psRegionSet(0,0,0,0)); // Covariance matrices
     1378    if (!images) {
     1379        psError(PS_ERR_IO, false, "Unable to read covariance matrices for chip %s, cell %s",
     1380                chipName, cellName);
     1381        return false;
     1382    }
     1383
     1384    psArray *readouts = cell->readouts; // Readouts of cell
     1385    if (images->n != readouts->n) {
     1386        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     1387                "Number of covariance matrices (%ld) doesn't match number of readouts (%ld)",
     1388                images->n, readouts->n);
     1389        psFree(images);
     1390        return false;
     1391    }
     1392
     1393    for (int i = 0; i < readouts->n; i++) {
     1394        pmReadout *ro = readouts->data[i]; // Readout of interest
     1395        psImage *image = images->data[i]; // Image of interest
     1396        if (ro->covariance) {
     1397            psWarning("Clobbering extant covariance matrix in chip %s, cell %s, readout %d",
     1398                      chipName, cellName, i);
     1399            psFree(ro->covariance);
     1400        }
     1401        ro->covariance = psKernelAllocFromImage(image, x0, y0);
     1402    }
     1403    psFree(images);
     1404
     1405    return true;
     1406}
     1407
     1408
     1409bool pmChipReadCovariance(pmChip *chip, psFits *fits)
     1410{
     1411    PS_ASSERT_PTR_NON_NULL(chip, false);
     1412    PS_ASSERT_FITS_NON_NULL(fits, false);
     1413
     1414    psArray *cells = chip->cells;       // Array of cells
     1415    for (int i = 0; i < cells->n; i++) {
     1416        pmCell *cell = cells->data[i];  // Cell of interest
     1417        if (!pmCellReadCovariance(cell, fits)) {
     1418            return false;
     1419        }
     1420    }
     1421
     1422    return true;
     1423}
     1424
     1425
     1426bool pmFPAReadCovariance(pmFPA *fpa, psFits *fits)
     1427{
     1428    PS_ASSERT_PTR_NON_NULL(fpa, false);
     1429    PS_ASSERT_FITS_NON_NULL(fits, false);
     1430
     1431    psArray *chips = fpa->chips;        // Array of chips
     1432    for (int i = 0; i < chips->n; i++) {
     1433        pmChip *chip = chips->data[i];  // Chip of interest
     1434        if (!pmChipReadCovariance(chip, fits)) {
     1435            return false;
     1436        }
     1437    }
     1438
     1439    return true;
     1440}

  • trunk/psModules/src/camera/pmFPARead.h

    r18163 r21363  
    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.16 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:24 $
    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
    … …  
    226226/// also read and included in the cell analysis metadata under "name.HEADER".
    227227int pmCellReadTable(pmCell *cell,       ///< Cell for which to read table
    228                     psFits *fits,       ///< FITS file from which the table
     228                    psFits *fits,       ///< FITS file from which to read the table
    229229                    const char *name    ///< Specifies the extension name, and target in the analysis metadata
    230230                   );
    … …  
    233233///
    234234/// Iterates over component cells, calling pmCellReadTable.
    235 int pmChipReadTable(pmChip *chip,       ///< Cell for which to read table
    236                     psFits *fits,       ///< FITS file from which the table
     235int pmChipReadTable(pmChip *chip,       ///< Chip for which to read table
     236                    psFits *fits,       ///< FITS file from which to read the table
    237237                    const char *name    ///< Specifies the extension name, and target in the analysis metadata
    238238                   );
    … …  
    241241///
    242242/// Iterates over component chips, calling pmChipReadTable.
    243 int pmFPAReadTable(pmFPA *fpa,          ///< Cell for which to read table
    244                    psFits *fits,        ///< FITS file from which the table
     243int pmFPAReadTable(pmFPA *fpa,          ///< FPA for which to read table
     244                   psFits *fits,        ///< FITS file from which to read the table
    245245                   const char *name     ///< Specifies the extension name, and target in the analysis metadata
    246246                  );
     247
     248/// Read covariance matrices for a cell
     249bool pmCellReadCovariance(pmCell *cell, ///< Cell for which to read covariance matrices
     250                          psFits *fits  ///< FITS file from which to read
     251    );
     252
     253/// Read covariance matrices for a chip
     254bool pmChipReadCovariance(pmChip *chip, ///< Chip for which to read covariance matrices
     255                          psFits *fits  ///< FITS file from which to read
     256    );
     257
     258/// Read covariance matrices for a cell
     259bool pmFPAReadCovariance(pmFPA *fpa,    ///< FPA for which to read covariance matrices
     260                         psFits *fits   ///< FITS file from which to read
     261    );
     262
     263
     264
    247265/// @}
    248266#endif

  • trunk/psModules/src/camera/pmFPAWrite.c

    r21279 r21363  
    99
    1010#include "pmConfig.h"
     11#include "pmConfigMask.h"
    1112#include "pmHDU.h"
    1213#include "pmFPA.h"
    … …  
    2930    FPA_WRITE_TYPE_IMAGE,               // Write image
    3031    FPA_WRITE_TYPE_MASK,                // Write mask
    31     FPA_WRITE_TYPE_WEIGHT               // Write weight map
     32    FPA_WRITE_TYPE_VARIANCE             // Write variance map
    3233} fpaWriteType;
    3334
    … …  
    4647    case FPA_WRITE_TYPE_MASK:
    4748        return &hdu->masks;
    48     case FPA_WRITE_TYPE_WEIGHT:
    49         return &hdu->weights;
     49    case FPA_WRITE_TYPE_VARIANCE:
     50        return &hdu->variances;
    5051    default:
    5152        psAbort("Unknown write type: %x\n", type);
    … …  
    6667    case FPA_WRITE_TYPE_MASK:
    6768        return pmHDUWriteMask(hdu, fits, config);
    68     case FPA_WRITE_TYPE_WEIGHT:
    69         return pmHDUWriteWeight(hdu, fits, config);
     69    case FPA_WRITE_TYPE_VARIANCE:
     70        return pmHDUWriteVariance(hdu, fits, config);
    7071    default:
    7172        psAbort("Unknown write type: %x\n", type);
    … …  
    7475}
    7576
    76 // Write a cell image/mask/weight
     77// Indicate whether a covariance matrix is defined
     78static bool readoutSearchCovariances(pmReadout *ro)
     79{
     80    return ro->covariance ? true : false;
     81}
     82
     83// Search for a covariance matrix
     84#define SEARCH_COVARIANCES(NAME, PARENT, CHILD, CHILDREN, TESTFUNC) \
     85static bool NAME(PARENT *parent) \
     86{ \
     87    if (!parent || !parent->CHILDREN) { \
     88        return false; \
     89    } \
     90    psArray *children = parent->CHILDREN; /* Array of children */ \
     91    for (int i = 0; i < children->n; i++) { \
     92        CHILD *child = children->data[i]; /* Child of interest */ \
     93        if (child && TESTFUNC(child)) { \
     94            return true; \
     95        } \
     96    } \
     97    return false; \
     98}
     99
     100SEARCH_COVARIANCES(cellSearchCovariances, pmCell, pmReadout, readouts, readoutSearchCovariances);
     101SEARCH_COVARIANCES(chipSearchCovariances, pmChip, pmCell,    cells,    cellSearchCovariances);
     102SEARCH_COVARIANCES(fpaSearchCovariances,  pmFPA,  pmChip,    chips,    chipSearchCovariances);
     103
     104// Some type-specific additions to the header
     105static bool writeUpdateHeader(pmFPA *fpa, // FPA of interest
     106                              pmChip *chip, // Chip of interest, or NULL
     107                              pmCell *cell, // Cell of interest, or NULL
     108                              fpaWriteType type, // Type to write
     109                              pmConfig *config // Configuration
     110                              )
     111{
     112    switch (type) {
     113      case FPA_WRITE_TYPE_MASK: {
     114          pmHDU *phu = pmHDUGetHighest(fpa, chip, cell); // Primary header
     115          if (!pmConfigMaskWriteHeader(config, phu->header)) {
     116              psError(PS_ERR_UNKNOWN, false, "Unable to set the mask names in the PHU header");
     117              return false;
     118          }
     119          break;
     120      }
     121      case FPA_WRITE_TYPE_VARIANCE: {
     122          bool covar = false;           // Are covariances present?
     123          if ((cell && cellSearchCovariances(cell)) ||
     124              (!cell && ((chip && chipSearchCovariances(chip)) ||
     125                         (!chip && fpa && fpaSearchCovariances(fpa))))) {
     126              covar = true;
     127          }
     128
     129          pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // Header being written
     130          psMetadataAddBool(hdu->header, PS_LIST_TAIL, PM_HDU_COVARIANCE_KEYWORD, PS_META_REPLACE,
     131                            "Is a covariance matrix present?", covar);
     132          break;
     133      }
     134      default:
     135        break;
     136    }
     137
     138    return true;
     139}
     140
     141
     142// Write a cell image/mask/variance
    77143static bool cellWrite(pmCell *cell,     // Cell to write
    78144                      psFits *fits,     // FITS file to which to write
    … …  
    97163    // generate the HDU, but only copies the structure.
    98164    if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) {
    99         psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.\n");
     165        psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.");
    100166        return false;
    101167    }
    … …  
    106172
    107173    if (writeBlank || writeImage) {
    108 
    109174        pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS |
    110175                                 PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE;
    111176        if (!pmConceptsWriteCell(cell, source, true, config)) {
    112             psError(PS_ERR_IO, false, "Unable to write concepts for cell.\n");
     177            psError(PS_ERR_IO, false, "Unable to write concepts for cell.");
     178            return false;
     179        }
     180        if (!writeUpdateHeader(NULL, NULL, cell, type, config)) {
     181            psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
    113182            return false;
    114183        }
    … …  
    123192}
    124193
    125 // Write a chip image/mask/weight
     194// Write a chip image/mask/variance
    126195static bool chipWrite(pmChip *chip,     // Chip to write
    127196                      psFits *fits,     // FITS file to which to write
    … …  
    162231                return false;
    163232            }
     233
     234            if (!writeUpdateHeader(NULL, chip, NULL, type, config)) {
     235                psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
     236                return false;
     237            }
     238
    164239            if (!appropriateWriteFunc(hdu, fits, config, type)) {
    165240                psError(PS_ERR_IO, false, "Unable to write HDU for chip.\n");
    … …  
    186261
    187262
    188 // Write an FPA image/mask/weight
     263// Write an FPA image/mask/variance
    189264static bool fpaWrite(pmFPA *fpa,        // FPA to write
    190265                     psFits *fits,      // FITS file to which to write
    … …  
    223298            if (!pmConceptsWriteFPA(fpa, source, true, config)) {
    224299                psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n");
     300                return false;
     301            }
     302            if (!writeUpdateHeader(fpa, NULL, NULL, type, config)) {
     303                psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
    225304                return false;
    226305            }
    … …  
    433512
    434513
    435 bool pmCellWriteWeight(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
     514bool pmCellWriteVariance(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
    436515{
    437516    PS_ASSERT_PTR_NON_NULL(cell, false);
    438517    PS_ASSERT_PTR_NON_NULL(fits, false);
    439     return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_WEIGHT);
    440 }
    441 
    442 bool pmChipWriteWeight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
     518    if (!cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_VARIANCE)) {
     519        return false;
     520    }
     521    if (!pmCellWriteCovariance(fits, cell)) {
     522        return false;
     523    }
     524    return true;
     525}
     526
     527bool pmChipWriteVariance(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
    443528{
    444529    PS_ASSERT_PTR_NON_NULL(chip, false);
    445530    PS_ASSERT_PTR_NON_NULL(fits, false);
    446     return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT);
    447 }
    448 
    449 bool pmFPAWriteWeight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
     531    if (!chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) {
     532        return false;
     533    }
     534    if (!pmChipWriteCovariance(fits, chip)) {
     535        return false;
     536    }
     537    return true;
     538}
     539
     540bool pmFPAWriteVariance(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
    450541{
    451542    PS_ASSERT_PTR_NON_NULL(fpa, false);
    452543    PS_ASSERT_PTR_NON_NULL(fits, false);
    453     return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT);
     544    if (!fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) {
     545        return false;
     546    }
     547    if (!pmFPAWriteCovariance(fits, fpa)) {
     548        return false;
     549    }
     550    return true;
    454551}
    455552
    … …  
    526623    return numWrite;
    527624}
     625
     626bool pmCellWriteCovariance(psFits *fits, const pmCell *cell)
     627{
     628    PS_ASSERT_PTR_NON_NULL(cell, false);
     629    PS_ASSERT_PTR_NON_NULL(fits, false);
     630
     631    int numCovar = 0;
     632    psArray *readouts = cell->readouts; // Array of readouts
     633    for (int i = 0; i < readouts->n; i++) {
     634        pmReadout *readout = readouts->data[i]; // The readout of interest
     635        if (readout && readout->covariance) {
     636            numCovar++;
     637        }
     638    }
     639    if (numCovar == 0) {
     640        return true;
     641    }
     642    if (numCovar != readouts->n) {
     643        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     644                "Number of covariances (%d) doesn't match number of readouts (%ld)",
     645                numCovar, readouts->n);
     646        return false;
     647    }
     648
     649    // Check size of covariances
     650    int xMinCovar = INT_MAX, xMaxCovar = INT_MIN, yMinCovar = INT_MAX, yMaxCovar = INT_MIN; // Size
     651    for (int i = 0; i < readouts->n; i++) {
     652        pmReadout *readout = readouts->data[i]; // The readout of interest
     653        psAssert(readout, "Should be defined.");
     654        psKernel *covar = readout->covariance; // Covariance matrix
     655        psAssert(covar, "Should be defined.");
     656        xMinCovar = PS_MIN(xMinCovar, covar->xMin);
     657        xMaxCovar = PS_MAX(xMaxCovar, covar->xMax);
     658        yMinCovar = PS_MIN(yMinCovar, covar->yMin);
     659        yMaxCovar = PS_MAX(yMaxCovar, covar->yMax);
     660    }
     661
     662    // Correct covariances to common size
     663    psArray *images = psArrayAlloc(numCovar); // Array of images
     664    for (int i = 0; i < readouts->n; i++) {
     665        pmReadout *readout = readouts->data[i]; // The readout of interest
     666        psAssert(readout, "Should be defined.");
     667        psKernel *covar = readout->covariance; // Covariance matrix
     668        psAssert(covar, "Should be defined.");
     669        int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax;// Size
     670        if (xMin == xMinCovar && xMax == xMaxCovar && yMin == yMinCovar && yMax == yMaxCovar) {
     671            images->data[i] = psMemIncrRefCounter(covar->image);
     672        } else {
     673            psImage *new = psImageAlloc(xMaxCovar - xMinCovar + 1, yMaxCovar - yMinCovar + 1, PS_TYPE_F32);
     674            psImageInit(new, 0);
     675            psImageOverlaySection(new, covar->image, xMinCovar - xMin, yMinCovar - yMin, "=");
     676            images->data[i] = new;
     677        }
     678    }
     679
     680    // Determine extension name
     681    const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip
     682    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
     683    psString extname = NULL;            // Extension name
     684    psStringAppend(&extname, "COVAR_%s_%s", chipName, cellName);
     685
     686    // Generate header
     687    pmHDU *hdu = pmHDUFromCell(cell);   // HDU for cell
     688    psMetadata *header = psMetadataCopy(NULL, hdu->header); // Header to write
     689    psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.X", PS_META_REPLACE,
     690                     "Centre of covariance matrix in x", -xMinCovar);
     691    psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.Y", PS_META_REPLACE,
     692                     "Centre of covariance matrix in y", -yMinCovar);
     693
     694    // Write images
     695    if (!psFitsWriteImageCube(fits, header, images, extname)) {
     696        psError(PS_ERR_UNKNOWN, false, "Unable to write covariances from chip %s, cell %s to extension %s",
     697                chipName, cellName, extname);
     698        psFree(extname);
     699        psFree(header);
     700        psFree(images);
     701        return 0;
     702    }
     703    psFree(extname);
     704    psFree(header);
     705    psFree(images);
     706
     707    return true;
     708}
     709
     710
     711bool pmChipWriteCovariance(psFits *fits, const pmChip *chip)
     712{
     713    PS_ASSERT_PTR_NON_NULL(chip, false);
     714    PS_ASSERT_PTR_NON_NULL(fits, false);
     715
     716    psArray *cells = chip->cells;       // Array of cells
     717    for (int i = 0; i < cells->n; i++) {
     718        pmCell *cell = cells->data[i];  // Cell of interest
     719        if (!pmCellWriteCovariance(fits, cell)) {
     720            return false;
     721        }
     722    }
     723
     724    return true;
     725}
     726
     727
     728bool pmFPAWriteCovariance(psFits *fits, const pmFPA *fpa)
     729{
     730    PS_ASSERT_PTR_NON_NULL(fpa, false);
     731    PS_ASSERT_PTR_NON_NULL(fits, false);
     732
     733    psArray *chips = fpa->chips;        // Array of chips
     734    for (int i = 0; i < chips->n; i++) {
     735        pmChip *chip = chips->data[i];  // Chip of interest
     736        if (!pmChipWriteCovariance(fits, chip)) {
     737            return false;
     738        }
     739    }
     740
     741    return true;
     742}

  • trunk/psModules/src/camera/pmFPAWrite.h

    r21279 r21363  
    44 * @author Paul Price, IfA
    55 *
    6  * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2009-02-04 02:39:36 $
     6 * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:24 $
     8 *
    89 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    910 */
    … …  
    2425                        psFits *fits,   ///<  FITS file to which to write
    2526                        int z           ///<  Image plane to write
    26                        );
     27    );
    2728
    2829/// Write a cell to a FITS file
    … …  
    3637                 pmConfig *config,      ///<  Configuration
    3738                 bool blank             ///<  Write a blank PHU?
    38                 );
     39    );
    3940
    4041/// Write a chip to a FITS file
    … …  
    4950                 bool blank,            ///<  Write a blank PHU?
    5051                 bool recurse           ///<  Recurse to lower levels?
    51                 );
     52    );
    5253
    5354/// Write an FPA to a FITS file
    … …  
    6263                bool blank,             ///<  Write a blank PHU?
    6364                bool recurse            ///<  Recurse to lower levels?
    64                );
     65    );
    6566
    6667/// Write a cell mask to a FITS file
    … …  
    7475                     pmConfig *config,  ///<  Configuration
    7576                     bool blank         ///<  Write a blank PHU?
    76                     );
     77    );
    7778
    7879/// Write a chip mask to a FITS file
    … …  
    8889                     bool blank,        ///<  Write a blank PHU?
    8990                     bool recurse       ///<  Recurse to lower levels?
    90                     );
     91    );
    9192
    9293/// Write an FPA mask to a FITS file
    … …  
    102103                    bool blank,         ///<  Write a blank PHU?
    103104                    bool recurse        ///<  Recurse to lower levels?
    104                    );
     105    );
    105106
    106 /// Write a cell weight to a FITS file
     107/// Write a cell variance to a FITS file
    107108///
    108 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     109/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    109110/// 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                       );
     111/// the HDU variance to the FITS file.  This function should be called at the beginning of the output cell
     112/// loop with blank=true in order to produce the correct file structure.
     113bool pmCellWriteVariance(pmCell *cell,    ///<  Cell to write
     114                         psFits *fits,    ///<  FITS file to which to write
     115                         pmConfig *config, ///<  Configuration
     116                         bool blank       ///<  Write a blank PHU?
     117    );
    117118
    118 /// Write a chip weight to a FITS file
     119/// Write a chip variance to a FITS file
    119120///
    120 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     121/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    121122/// 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
     123/// the HDU variance to the FITS file, optionally recursing to lower levels.  This function should be called
     124/// at the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct
    124125/// 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                       );
     126bool pmChipWriteVariance(pmChip *chip,    ///<  Chip to write
     127                         psFits *fits,    ///<  FITS file to which to write
     128                         pmConfig *config, ///<  Configuration
     129                         bool blank,      ///<  Write a blank PHU?
     130                         bool recurse     ///<  Recurse to lower levels?
     131    );
    131132
    132 /// Write an FPA weight to a FITS file
     133/// Write an FPA variance to a FITS file
    133134///
    134 /// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
     135/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
    135136/// 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
     137/// the HDU variance to the FITS file, optionally recursing to lower levels.  This function should be called
     138/// at the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct
    138139/// 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                      );
     140bool pmFPAWriteVariance(pmFPA *fpa,       ///<  FPA to write
     141                        psFits *fits,     ///<  FITS file to which to write
     142                        pmConfig *config, ///<  Configuration
     143                        bool blank,       ///<  Write a blank PHU?
     144                        bool recurse      ///<  Recurse to lower levels?
     145    );
    145146
    146147
    … …  
    153154                     const pmCell *cell, ///< Cell containing FITS table in the analysis metadata
    154155                     const char *name   ///< Name for the table data, and the extension name
    155                     );
     156    );
    156157
    157158int pmChipWriteTable(psFits *fits,      ///< FITS file to which to write
    158159                     const pmChip *chip, ///< Chip containing cells with tables to write
    159160                     const char *name   ///< Name for the table data, and the extension name
    160                     );
     161    );
    161162
    162163int pmFPAWriteTable(psFits *fits,       ///< FITS file to which to write
    163164                    const pmFPA *fpa,   ///< FPA containing cells with tables to write
    164165                    const char *name    ///< Name for the table data, and the extension name
    165                    );
     166    );
     167
     168/// Write covariance matrix to a FITS file
     169///
     170/// The covariance matrices for a cell are written to an independent extension, named after the chip and cell
     171/// name.
     172bool pmCellWriteCovariance(psFits *fits,///< FITS file to which to write
     173                           const pmCell *cell ///< Cell for which to write covariance
     174    );
     175
     176bool pmChipWriteCovariance(psFits *fits,///< FITS file to which to write
     177                           const pmChip *chip ///< Chip for which to write covariance
     178    );
     179
     180bool pmFPAWriteCovariance(psFits *fits,///< FITS file to which to write
     181                          const pmFPA *fpa ///< FPA for which to write covariance
     182    );
    166183
    167184// Update the header before writing to be consistent

  • trunk/psModules/src/camera/pmFPAfile.c

    r21279 r21363  
    476476        return PM_FPA_FILE_MASK;
    477477    }
    478     if (!strcasecmp(type, "WEIGHT"))     {
    479         return PM_FPA_FILE_WEIGHT;
     478    if (!strcasecmp(type, "VARIANCE"))     {
     479        return PM_FPA_FILE_VARIANCE;
    480480    }
    481481    if (!strcasecmp(type, "FRINGE")) {
    … …  
    530530      case PM_FPA_FILE_MASK:
    531531        return ("MASK");
    532       case PM_FPA_FILE_WEIGHT:
    533         return ("WEIGHT");
     532      case PM_FPA_FILE_VARIANCE:
     533        return ("VARIANCE");
    534534      case PM_FPA_FILE_FRINGE:
    535535        return ("FRINGE");

  • trunk/psModules/src/camera/pmFPAfile.h

    r21279 r21363  
    44 * @author EAM, IfA
    55 *
    6  * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2009-02-04 02:39:36 $
     6 * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:24 $
     8 *
    89 * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
    910 */
    … …  
    3738    PM_FPA_FILE_IMAGE,
    3839    PM_FPA_FILE_MASK,
    39     PM_FPA_FILE_WEIGHT,
     40    PM_FPA_FILE_VARIANCE,
    4041    PM_FPA_FILE_FRINGE,
    4142    PM_FPA_FILE_DARK,

  • trunk/psModules/src/camera/pmFPAfileFitsIO.c

    r21279 r21363  
    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

  • trunk/psModules/src/camera/pmFPAfileFitsIO.h

    r18601 r21363  
    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.16 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-02-06 02:31:24 $
    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

  • trunk/psModules/src/camera/pmFPAfileIO.c

    r20937 r21363  
    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:
    … …  
    771771          }
    772772
    773           // XXX if we have a mask file, then we need to read the mask bit names
    774           // defined for this file
    775           if (file->type == PM_FPA_FILE_MASK) {
    776             if (!pmConfigMaskReadHeader (config, phu)) {
    777                 psError(PS_ERR_IO, false, "error in mask bits");
    778                 return false;
    779             }
    780           }
    781 
    782773          // determine the current format from the header
    783774          // determine camera if not specified already
    784775          // XXX can I actually reach this with camera not specified??
    785           psMetadata *camera = NULL;
    786           psString formatName = NULL;
    787           psString cameraName = NULL;
    788           file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true);
     776          psMetadata *camera = NULL;
     777          psString formatName = NULL;
     778          psString cameraName = NULL;
     779          file->format = pmConfigCameraFormatFromHeader(&camera, &cameraName, &formatName, config, phu, true);
    789780          if (!file->format) {
    790781            psError(PS_ERR_IO, false, "Failed to read CCD format from %s\n", file->filename);
    … …  
    794785          psFree(phu);
    795786
    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;
     787          pmFPA *newFPA = pmFPAConstruct (camera, formatName);
     788          if (!newFPA) {
     789              psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName);
     790              psFree(camera);
     791              psFree(formatName);
     792              return NULL;
     793          }
     794          psFree(camera);
     795          psFree(formatName);
     796          psFree(cameraName);
     797
     798          // XXX this is really dangerous...
     799          psFree (file->fpa);
     800          file->fpa = newFPA;
    810801        }
    811802
    … …  
    938929      case PM_FPA_FILE_IMAGE:
    939930      case PM_FPA_FILE_MASK:
    940       case PM_FPA_FILE_WEIGHT:
     931      case PM_FPA_FILE_VARIANCE:
    941932      case PM_FPA_FILE_DARK:
    942933      case PM_FPA_FILE_FRINGE:

  • trunk/psModules/src/camera/pmHDU.c

    r21313 r21363  
    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);
     243    return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits);
    244244}
    245245

  • trunk/psModules/src/camera/pmHDU.h

    r21279 r21363  
    44 * @author Paul Price, IfA
    55 *
    6  * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2009-02-04 02:39:36 $
     6 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:24 $
     8 *
    89 * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
    910 */
    … …  
    1819/// @{
    1920
     21#define PM_HDU_COVARIANCE_KEYWORD "PS_COVAR" // FITS keyword to indicate presence of a covariance matrix
     22
    2023/// An instance of the FITS Header Data Unit
    2124///
    22 /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weight data, but these are
     25/// Of course, it is not an exact replica of a FITS HDU --- they have no mask and variance data, but these are
    2326/// stored here for convenience --- it keeps all the relevant data about the image in one place.
    24 typedef struct
    25 {
     27typedef struct {
    2628    psString extname;                   ///< The extension name
    2729    bool blankPHU;                      ///< Is this a blank FITS Primary Header Unit, i.e., no data?
    2830    psMetadata *format;                 ///< The camera format
    2931    psMetadata *header;                 ///< The FITS header, or NULL if primary for FITS; or section info
    30     psArray *images;                    ///< The pixel data
    31     psArray *weights;                   ///< The pixel data
    32     psArray *masks;                     ///< The pixel data
    33 }
    34 pmHDU;
     32    psArray *images;                    ///< Pixel data
     33    psArray *variances;                 ///< Variance in the pixel data, or NULL
     34    psArray *masks;                     ///< Mask for the pixel data, or NULL
     35} pmHDU;
    3536
    3637
    … …  
    6061                  );
    6162
    62 /// Read the HDU header and weight map
     63/// Read the HDU header and variance map
    6364///
    6465/// Moves to the appropriate extension
    65 bool pmHDUReadWeight(pmHDU *hdu,        ///< HDU to read
    66                      psFits *fits       ///< FITS file to read from
     66bool pmHDUReadVariance(pmHDU *hdu,        ///< HDU to read
     67                       psFits *fits       ///< FITS file to read from
    6768    );
    6869
    … …  
    7980    );
    8081
    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
     82/// Write the HDU header and variance map
     83bool pmHDUWriteVariance(pmHDU *hdu,       ///< HDU to write
     84                        psFits *fits,     ///< FITS file to write to
     85                        const pmConfig *config  ///< Configuration
    8586    );
    8687

  • trunk/psModules/src/camera/pmHDUGenerate.c

    r18227 r21363  
    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    }
    … …  
    448448
    449449            psArray *readouts = cell->readouts; // Array of readouts
     450
    450451            psArray *hduImages = hdu->images; // Array of images in the HDU
    451452            psArray *hduMasks = hdu->masks; // Array of masks in the HDU
    452             psArray *hduWeights = hdu->weights; // Array of weights in the HDU
     453            psArray *hduVariances = hdu->variances; // Array of variances in the HDU
    453454            for (int i = 0; i < readouts->n; i++) {
    454455                pmReadout *readout = readouts->data[i]; // The readout of interest
    … …  
    467468                    readout->mask = new;
    468469                }
    469                 if (readout->weight) {
    470                     psImage *new = pasteImage(hduWeights->data[i], readout->weight, trimsec);
    471                     psFree(readout->weight);
    472                     readout->weight = new;
     470                if (readout->variance) {
     471                    psImage *new = pasteImage(hduVariances->data[i], readout->variance, trimsec);
     472                    psFree(readout->variance);
     473                    readout->variance = new;
    473474                }
    474475
    … …  
    504505
    505506// Return the level that an extension applies to
    506 static pmFPALevel extensionLevel(pmHDU *hdu // HDU to check
     507static pmFPALevel extensionLevel(const pmHDU *hdu // HDU to check
    507508                                )
    508509{
    … …  
    512513    }
    513514    bool mdok = true;                   // Status of MD lookup
    514     psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File information for camera format
     515    psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File info for camera format
    515516    if (!mdok || !file) {
    516517        psError(PS_ERR_UNEXPECTED_NULL, true, "Can't file FILE information for camera format "
    … …  
    581582        return false;
    582583    }
    583     if (hdu->images && hdu->masks && hdu->weights) {
     584    if (hdu->images && hdu->masks && hdu->variances) {
    584585        // It's already here!
    585586        return true;
    … …  
    631632        return generateForCells(chip);
    632633    }
    633     if (hdu->images && hdu->masks && hdu->weights) {
     634    if (hdu->images && hdu->masks && hdu->variances) {
    634635        // It's already here!
    635636        return true;
    … …  
    680681        return generateForChips(fpa);
    681682    }
    682     if (hdu->images && hdu->masks && hdu->weights) {
     683    if (hdu->images && hdu->masks && hdu->variances) {
    683684        // It's already here!
    684685        return true;

  • trunk/psModules/src/camera/pmHDUUtils.c

    r18554 r21363  
    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

  • trunk/psModules/src/camera/pmReadoutStack.c

    r21183 r21363  
    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    }

  • trunk/psModules/src/config/pmConfigMask.c

    r21183 r21363  
    1010
    1111static pmConfigMaskInfo masks[] = {
    12     { "DETECTOR",  NULL,       0x00, true  },   // Something is wrong with the detector
    13     { "DARK",      "DETECTOR", 0x00, true  },   // Pixel doesn't dark-subtract properly
    14     { "FLAT",      "DETECTOR", 0x01, true  },   // Pixel doesn't flat-field properly
    15     { "BLANK",     "DETECTOR", 0x01, true  },   // Pixel doesn't contain valid data
    16     { "RANGE",     NULL,       0x00, true  },   // Pixel is out-of-range of linearity
    17     { "SAT",       "RANGE",    0x01, true  },   // Pixel is saturated
    18     { "BAD",       "RANGE",    0x01, true  },   // Pixel is low
    19     { "BAD.WARP",  NULL,       0x01, true  },   // Pixel is bad after convolution with a bad pixel
    20     { "CR",        NULL,       0x00, true  },   // Pixel contains a cosmic ray
    21     { "GHOST",     NULL,       0x00, true  },   // Pixel contains an optical ghost
    22     { "POOR.WARP", NULL,       0x00, false },   // Pixel is poor after convolution with a bad pixel
     12    { "DETECTOR",  NULL,       0x00, true  },   // Something is wrong with the detector
     13    { "DARK",      "DETECTOR", 0x00, true  },   // Pixel doesn't dark-subtract properly
     14    { "FLAT",      "DETECTOR", 0x01, true  },   // Pixel doesn't flat-field properly
     15    { "BLANK",     "DETECTOR", 0x01, true  },   // Pixel doesn't contain valid data
     16    { "RANGE",     NULL,       0x00, true  },   // Pixel is out-of-range of linearity
     17    { "SAT",       "RANGE",    0x01, true  },   // Pixel is saturated
     18    { "BAD",       "RANGE",    0x01, true  },   // Pixel is low
     19    { "BAD.WARP",  NULL,       0x01, true  },   // Pixel is bad after convolution with a bad pixel
     20    { "CR",        NULL,       0x00, true  },   // Pixel contains a cosmic ray
     21    { "GHOST",     NULL,       0x00, true  },   // Pixel contains an optical ghost
     22    { "POOR.WARP", NULL,       0x00, false },   // Pixel is poor after convolution with a bad pixel
    2323    // "LOW"  Pixel is low
    2424    // "CONV" Pixel is bad after convolution with a bad pixel
    … …  
    3131// bits in 8-bits of space).
    3232
    33 // XXX this file does not have psError vs psWarning worked out correctly.  some of the 
     33// XXX this file does not have psError vs psWarning worked out correctly.  some of the
    3434// failure modes should result in errors, not just warnings.
    3535
    … …  
    4040
    4141bool pmConfigMaskSetInMetadata(psImageMaskType *outMaskValue, // Value of MASK.VALUE, returned
    42                                psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned
    43                                psMetadata *source  // Source of mask bits
     42                               psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned
     43                               psMetadata *source  // Source of mask bits
    4444    )
    4545{
    4646    PS_ASSERT_METADATA_NON_NULL(source, false);
    47    
     47
    4848    // Ensure all the bad mask names exist, and set the value to catch all bad pixels
    4949    psImageMaskType maskValue = 0;           // Value to mask to catch all the bad pixels
    50     psImageMaskType allMasks = 0;            // Value to mask to catch all masked bits (to set MARK)
     50    psImageMaskType allMasks = 0;            // Value to mask to catch all masked bits (to set MARK)
    5151
    5252    int nMasks = sizeof (masks) / sizeof (pmConfigMaskInfo);
    … …  
    5555        bool mdok;                      // Status of MD lookup
    5656        psImageMaskType value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].badMaskName); // Value of mask
    57         if (!mdok) {
    58             psWarning ("problem with mask value %s\n", masks[i].badMaskName);
    59         }
     57        if (!mdok) {
     58            psWarning ("problem with mask value %s\n", masks[i].badMaskName);
     59        }
    6060
    6161        if (!value) {
    62             if (masks[i].fallbackName) {
    63                 value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].fallbackName);
    64             }
    65             if (!value) {
    66                 value = masks[i].defaultMaskValue;
    67             }
    68             psMetadataAddImageMask(source, PS_LIST_TAIL, masks[i].badMaskName, PS_META_REPLACE, NULL, value);
    69         }
    70         if (masks[i].isBad) {
    71             maskValue |= value;
    72         }
    73         allMasks |= value;
     62            if (masks[i].fallbackName) {
     63                value = psMetadataLookupImageMaskFromGeneric(&mdok, source, masks[i].fallbackName);
     64            }
     65            if (!value) {
     66                value = masks[i].defaultMaskValue;
     67            }
     68            psMetadataAddImageMask(source, PS_LIST_TAIL, masks[i].badMaskName, PS_META_REPLACE, NULL, value);
     69        }
     70        if (masks[i].isBad) {
     71            maskValue |= value;
     72        }
     73        allMasks |= value;
    7474    }
    7575
    … …  
    107107// Get a mask value by name(s)
    108108psImageMaskType pmConfigMaskGetFromMetadata(psMetadata *source, // Source of masks
    109                                             const char *masks // Mask values to get
     109                                            const char *masks // Mask values to get
    110110    )
    111111{
    … …  
    139139}
    140140
    141 // lookup an image mask value by name from a psMetadata, without requiring the entry to 
     141// lookup an image mask value by name from a psMetadata, without requiring the entry to
    142142// be of type psImageMaskType, but verifying that it will fit in psImageMaskType
    143 psImageMaskType psMetadataLookupImageMaskFromGeneric (bool *status, const psMetadata *md, const char *name) {
    144 
     143psImageMaskType psMetadataLookupImageMaskFromGeneric(bool *status, const psMetadata *md, const char *name)
     144{
     145    if (!md) {
     146        psError(PS_ERR_UNEXPECTED_NULL, true, "Metadata is NULL.");
     147        if (status) {
     148            *status = false;
     149        }
     150        return 0;
     151    }
     152    if (!name) {
     153        psError(PS_ERR_UNEXPECTED_NULL, true, "Keyword is NULL.");
     154        if (status) {
     155            *status = false;
     156        }
     157        return 0;
     158    }
    145159    *status = true;
    146160
    147         // select the mask bit name from the header
    148         psMetadataItem *item = psMetadataLookup (md, name);
    149         if (!item) {
    150             psWarning("Unable to find header keyword %s when parsing mask", name);
    151             *status = false;
    152             return 0;
    153         }
    154 
    155         // the value may be any of the U8, U16, U32, U64 types : accept the value regardless of type size
    156         psU64 fullValue = 0;
    157         switch (item->type) {
    158           case PS_DATA_U8:
    159             fullValue = item->data.U8;
    160             break;
    161           case PS_DATA_U16:
    162             fullValue = item->data.U16;
    163             break;
    164           case PS_DATA_U32:
    165             fullValue = item->data.U32;
    166             break;
    167           case PS_DATA_U64:
    168             fullValue = item->data.U64;
    169             break;
    170           case PS_DATA_S8:
    171             fullValue = item->data.S8;
    172             break;
    173           case PS_DATA_S16:
    174             fullValue = item->data.S16;
    175             break;
    176           case PS_DATA_S32:
    177             fullValue = item->data.S32;
    178             break;
    179           case PS_DATA_S64:
    180             fullValue = item->data.S64;
    181             break;
    182           default:
    183             psWarning("Mask entry %s in metadata is not of a mask type", name);
    184             *status = false;
    185             return 0;
    186         }
    187 
    188         // will the incoming value fit within the current image mask type?
    189         if (fullValue > PS_MAX_IMAGE_MASK_TYPE) {
    190             psWarning("Mask entry %s in metadata is larger than allowed by the psImageMaskType", name);
    191             *status = false;
    192             return 0;
    193         }
    194         psImageMaskType value = fullValue;
    195         // XXX validate that value is a 2^n value?
    196 
    197         return value;
     161    // select the mask bit name from the header
     162    psMetadataItem *item = psMetadataLookup(md, name);
     163    if (!item) {
     164        if (status) {
     165            *status = false;
     166        } else {
     167            psError(PS_ERR_BAD_PARAMETER_VALUE, "Unable to find keyword %s when parsing mask", name);
     168        }
     169        return 0;
     170    }
     171
     172    // the value may be any of the U8, U16, U32, U64 types : accept the value regardless of type size
     173    psU64 fullValue = 0;
     174    switch (item->type) {
     175      case PS_DATA_U8:
     176        fullValue = item->data.U8;
     177        break;
     178      case PS_DATA_U16:
     179        fullValue = item->data.U16;
     180        break;
     181      case PS_DATA_U32:
     182        fullValue = item->data.U32;
     183        break;
     184      case PS_DATA_U64:
     185        fullValue = item->data.U64;
     186        break;
     187      case PS_DATA_S8:
     188        fullValue = item->data.S8;
     189        break;
     190      case PS_DATA_S16:
     191        fullValue = item->data.S16;
     192        break;
     193      case PS_DATA_S32:
     194        fullValue = item->data.S32;
     195        break;
     196      case PS_DATA_S64:
     197        fullValue = item->data.S64;
     198        break;
     199      default:
     200        if (status) {
     201            *status = false;
     202        } else {
     203            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     204                    "Mask entry %s in metadata is not of a mask type", name);
     205        }
     206        return 0;
     207    }
     208
     209    // will the incoming value fit within the current image mask type?
     210    if (fullValue > PS_MAX_IMAGE_MASK_TYPE) {
     211        if (status) {
     212            *status = false;
     213        } else {
     214            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     215                    "Mask entry %s in metadata is larger than allowed by the psImageMaskType", name);
     216        }
     217        return 0;
     218    }
     219    psImageMaskType value = fullValue;
     220    // XXX validate that value is a 2^n value?
     221
     222    return value;
    198223}
    199224
    200225// Remove from the header keywords starting with the provided string
    201226int pmConfigMaskRemoveHeaderKeywords(psMetadata *header, // Header from which to remove keywords
    202                                      const char *start // Remove keywords that start with this string
     227                                     const char *start // Remove keywords that start with this string
    203228    )
    204229{
    … …  
    229254        return 0;
    230255    }
    231    
     256
    232257    psImageMaskType mask = pmConfigMaskGetFromMetadata (recipe, masks);
    233258    return mask;
    … …  
    274299    }
    275300
    276     // How many mask values do we need to read?  We raise an error if this is not found, 
     301    // How many mask values do we need to read?  We raise an error if this is not found,
    277302    // unless the MASK.FORCE is set to true in the camera config
    278303    int nMask = psMetadataLookupS32(&status, header, "MSKNUM");
    … …  
    300325        }
    301326
    302         psImageMaskType headerValue = psMetadataLookupImageMaskFromGeneric (&status, header, valuekey);
    303         if (!status) {
     327        psImageMaskType headerValue = psMetadataLookupImageMaskFromGeneric (&status, header, valuekey);
     328        if (!status) {
    304329            psWarning("Failed to get mask value %s from header, skipping", valuekey);
    305             continue;
    306         }           
    307 
    308         // since we may read multiple mask files, we need to warn (or error?) if any of the
    309         // header mask values conflict with other header mask values; However, the original
    310         // mask values from the recipe do not need to match the header values.
    311 
    312         // when we add a header mask value, we will also add the NAME.ALREADY entry; check for
    313         // the NAME.ALREADY entry to see if we have previously added this mask value from a
    314         // header.
     330            continue;
     331        }
     332
     333        // since we may read multiple mask files, we need to warn (or error?) if any of the
     334        // header mask values conflict with other header mask values; However, the original
     335        // mask values from the recipe do not need to match the header values.
     336
     337        // when we add a header mask value, we will also add the NAME.ALREADY entry; check for
     338        // the NAME.ALREADY entry to see if we have previously added this mask value from a
     339        // header.
    315340
    316341        psString nameAlready = NULL;    // Name of key with ".ALREADY" added
    … …  
    318343        bool already = psMetadataLookupBool(&status, recipe, nameAlready); // Already read this one?
    319344
    320         bool inRecipe = false;
    321         psImageMaskType recipeValue = psMetadataLookupImageMaskFromGeneric (&inRecipe, recipe, name);
    322         if (!inRecipe) {
     345        bool inRecipe = false;
     346        psImageMaskType recipeValue = psMetadataLookupImageMaskFromGeneric (&inRecipe, recipe, name);
     347        if (!inRecipe) {
    323348            psWarning("Mask value %s is not defined in the recipe", name);
    324         }           
     349        }
    325350
    326351        if (already) {
    327             assert (inRecipe); // XXX makes no sense for NAME.ALREADY to be in without NAME
     352            assert (inRecipe); // XXX makes no sense for NAME.ALREADY to be in without NAME
    328353            if (recipeValue != headerValue) {
    329354                psWarning("New mask header value does not match previously loaded entry: %x vs %x", headerValue, recipeValue);
    330                 psMetadataAddImageMask(recipe, PS_LIST_TAIL, name, PS_META_REPLACE, "Bitmask bit value", headerValue);
    331                 // XXX alternatively, error here
     355                psMetadataAddImageMask(recipe, PS_LIST_TAIL, name, PS_META_REPLACE, "Bitmask bit value", headerValue);
     356                // XXX alternatively, error here
    332357            }
    333358        } else {
    … …  
    369394    while ((item = psMetadataGetAndIncrement(iter))) {
    370395
    371         // XXX this would give a false positive for mask which include '.ALREADY' in their names
    372         char *ptr = strstr (item->name, ".ALREADY");
    373         if (ptr) continue;
    374 
    375         psU64 fullValue = 0;
     396        // XXX this would give a false positive for mask which include '.ALREADY' in their names
     397        char *ptr = strstr (item->name, ".ALREADY");
     398        if (ptr) continue;
     399
     400        psU64 fullValue = 0;
    376401        switch (item->type) {
    377           case PS_DATA_U8:
    378             fullValue = item->data.U8;
    379             break;
    380           case PS_DATA_U16:
    381             fullValue = item->data.U16;
    382             break;
    383           case PS_DATA_U32:
    384             fullValue = item->data.U32;
    385             break;
    386           case PS_DATA_U64:
    387             fullValue = item->data.U64;
    388             break;
    389           default:
     402          case PS_DATA_U8:
     403            fullValue = item->data.U8;
     404            break;
     405          case PS_DATA_U16:
     406            fullValue = item->data.U16;
     407            break;
     408          case PS_DATA_U32:
     409            fullValue = item->data.U32;
     410            break;
     411          case PS_DATA_U64:
     412            fullValue = item->data.U64;
     413            break;
     414          default:
    390415            psWarning("mask recipe entry %s is not a bit value\n", item->name);
    391416            continue;
    392417        }
    393         assert (fullValue <= PS_MAX_IMAGE_MASK_TYPE); // this should have been asserted on read...
     418        assert (fullValue <= PS_MAX_IMAGE_MASK_TYPE); // this should have been asserted on read...
    394419
    395420        snprintf(namekey,  64, "MSKNAM%02d", nMask);

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r21183 r21363  
    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

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r21183 r21363  
    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);

  • trunk/psModules/src/imcombine/pmReadoutCombine.c

    r21183 r21363  
    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);

  • trunk/psModules/src/imcombine/pmReadoutCombine.h

    r21183 r21363  
    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.16 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-02-06 02:31:25 $
    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

  • trunk/psModules/src/imcombine/pmStack.c

    r21183 r21363  
    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.47 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2009-02-06 02:31:25 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    3030#define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
    3131#define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
    32 //#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
     32#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
     33//#define ADD_VARIANCE                  // Allow additional variance (besides variance factor)?
    3334#define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
    3435
    … …  
    263264
    264265        psImage *image = data->readout->image; // Image of interest
    265         psImage *variance = data->readout->weight; // Variance ("weight") map of interest
     266        psImage *variance = data->readout->variance; // Variance map of interest
    266267        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
    267268        if (variance) {
    … …  
    442443    *numCols = data->readout->image->numCols;
    443444    *numRows = data->readout->image->numRows;
    444     if (data->readout->weight) {
     445    if (data->readout->variance) {
    445446        *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);
     447        PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false);
     448        PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false);
     449        PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
    449450    }
    450451    *haveRejects = (data->reject != NULL);
    … …  
    473474        PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false);
    474475        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);
     476            PS_ASSERT_IMAGE_NON_NULL(data->readout->variance, false);
     477            PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->variance, false);
     478            PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
    478479        }
    479480    }
    … …  
    606607        weights->data.F32[i] = data->weight;
    607608        stack->data[i] = psMemIncrRefCounter(data->readout);
    608         // Variance factor
    609         float vf = psMetadataLookupF32(NULL, data->readout->parent->concepts, "CELL.VARFACTOR"); // Var factor
    610         if (!isfinite(vf)) {
    611             psWarning("Non-finite CELL.VARFACTOR for image %d --- setting to unity.", i);
    612             vf = 1.0;
    613         }
    614         varFactors->data.F32[i] = vf;
     609        varFactors->data.F32[i] = data->readout->covariance ?
     610            psImageCovarianceFactor(data->readout->covariance) : 1.0;
     611#if 0
    615612        if (isfinite(data->addVariance)) {
    616613            varFactors->data.F32[i] *= data->addVariance;
    617614        }
     615#endif
    618616        if (!haveRejects && !data->inspect) {
    619617            data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
    … …  
    649647        psImage *combinedImage = combined->image; // Combined image
    650648        psImage *combinedMask = combined->mask; // Combined mask
    651         psImage *combinedVariance = combined->weight; // Combined variance map
     649        psImage *combinedVariance = combined->variance; // Combined variance map
    652650
    653651        psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
    … …  
    702700        }
    703701
    704         psImage *combinedVariance = combined->weight; // Combined variance map
     702        psImage *combinedVariance = combined->variance; // Combined variance map
    705703        if (haveVariances && !combinedVariance) {
    706             combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    707             combinedVariance = combined->weight;
     704            combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     705            combinedVariance = combined->variance;
    708706        }
    709707

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r21351 r21363  
    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);
    … …  
    915915}
    916916
    917 psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
     917psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
    918918{
    919919    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
    … …  
    922922    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
    923923
    924     // Precalulate polynomial values
    925     psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y);
    926 
    927     // The appropriate kernel
    928     psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
    929 
     924    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial
     925    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel
    930926    psFree(polyValues);
    931927
     928    return kernel;
     929}
     930
     931psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
     932{
     933    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
     934    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
     935    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
     936    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
     937
     938    psKernel *kernel = pmSubtractionKernel(kernels, x, y, wantDual); // Convolution kernel
    932939    psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel
    933940    psFree(kernel);
    … …  
    989996
    990997
    991 // XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
     998// XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data
    992999static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
    9931000                                     int x0, int y0, // Offsets for image
    … …  
    10101017
    10111018    psKernel *kernelImage = NULL;       // Kernel for the images
    1012     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
     1019    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
    10131020
    10141021    // Only generate polynomial values every kernel footprint, since we have already assumed
    … …  
    10201027
    10211028    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,
     1029        convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance,
     1030                       ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,
    10241031                       maskBad, maskPoor, poorFrac, useFFT, false);
    10251032    }
    10261033    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,
     1034        convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance,
     1035                       ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region,
    10291036                       maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
    10301037    }
    10311038
    10321039    psFree(kernelImage);
    1033     psFree(kernelWeight);
     1040    psFree(kernelVariance);
    10341041    psFree(polyValues);
    10351042
    … …  
    11481155            // XXX }
    11491156        }
    1150         if (ro1->weight) {
    1151             if (!out1->weight) {
    1152                 out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     1157        if (ro1->variance) {
     1158            if (!out1->variance) {
     1159                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    11531160                // XXX if (threaded) {
    1154                 // XXX     psMutexInit(out1->weight);
     1161                // XXX     psMutexInit(out1->variance);
    11551162                // XXX }
    11561163            }
    1157             psImageInit(out1->weight, 0.0);
     1164            psImageInit(out1->variance, 0.0);
    11581165        }
    11591166    }
    … …  
    11651172            // XXX }
    11661173        }
    1167         if (ro2->weight) {
    1168             if (!out2->weight) {
    1169                 out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     1174        if (ro2->variance) {
     1175            if (!out2->variance) {
     1176                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    11701177                // XXX if (threaded) {
    1171                 // XXX     psMutexInit(out2->weight);
     1178                // XXX     psMutexInit(out2->variance);
    11721179                // XXX }
    11731180            }
    1174             psImageInit(out2->weight, 0.0);
     1181            psImageInit(out2->variance, 0.0);
    11751182        }
    11761183    }
    … …  
    12321239    psImage *polyValues = NULL;         // Pre-calculated polynomial values
    12331240    psKernel *kernelImage = NULL;       // Kernel for the images
    1234     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
     1241    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
    12351242#endif
    12361243
    … …  
    12991306                psFree(job);
    13001307            } else {
    1301                 subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
    1302                                          sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion,
    1303                                          kernels, doBG, useFFT);
     1308                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
     1309                                         subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,
     1310                                         useFFT);
    13041311            }
    13051312            psFree(subRegion);
    … …  
    13341341        // XXX }
    13351342    }
    1336 
    13371343    psImageConvolveSetThreads(oldThreads);
    13381344
    13391345    psFree(sys1);
    13401346    psFree(sys2);
     1347
     1348    // Calculate covariances
     1349    // This can take a while, so we only do it for a single instance
     1350    // XXX psImageCovarianceCalculate could be multithreaded
     1351    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     1352        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     1353        out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance);
     1354        psFree(kernel);
     1355    }
     1356    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     1357        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0,
     1358                                               kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel
     1359        out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance);
     1360        psFree(kernel);
     1361    }
    13411362
    13421363    // Copy anything that wasn't convolved
    … …  
    13451366        if (out2) {
    13461367            out2->image = psMemIncrRefCounter(ro2->image);
    1347             out2->weight = psMemIncrRefCounter(ro2->weight);
     1368            out2->variance = psMemIncrRefCounter(ro2->variance);
    13481369            out2->mask = psMemIncrRefCounter(ro2->mask);
     1370            out2->covariance = psMemIncrRefCounter(ro2->covariance);
    13491371        }
    13501372        break;
    … …  
    13521374        if (out1) {
    13531375            out1->image = psMemIncrRefCounter(ro1->image);
    1354             out1->weight = psMemIncrRefCounter(ro1->weight);
     1376            out1->variance = psMemIncrRefCounter(ro1->variance);
    13551377            out1->mask = psMemIncrRefCounter(ro1->mask);
     1378            out1->covariance = psMemIncrRefCounter(ro1->covariance);
    13561379        }
    13571380        break;

  • trunk/psModules/src/imcombine/pmSubtraction.h

    r21183 r21363  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2009-01-27 06:39:38 $
     8 * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2009-02-06 02:31:25 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    7070                              float sigmaRej, ///< Number of RMS deviations above zero at which to reject
    7171                              int footprint ///< Half-size of stamp
     72    );
     73
     74/// Generate the convolution kernel
     75psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, ///< Kernel parameters
     76                              float x, float y, ///< Normalised position [-1,1] for which to generate image
     77                              bool wantDual ///< Calculate for the dual kernel?
    7278    );
    7379

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r20561 r21363  
    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        }

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r21247 r21363  
    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;

  • trunk/psModules/src/imcombine/pmSubtractionParams.c

    r18287 r21363  
    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;

  • trunk/psModules/src/imcombine/pmSubtractionStamps.c

    r21183 r21363  
    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

  • trunk/psModules/src/imcombine/pmSubtractionStamps.h

    r20465 r21363  
    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    );

  • trunk/psModules/src/imcombine/pmSubtractionThreads.c

    r21351 r21363  
    2727    // XXX     psMutexInit(ro->image);
    2828    // XXX }
    29     // XXX if (ro->weight) {
    30     // XXX     psMutexInit(ro->weight);
     29    // XXX if (ro->variance) {
     30    // XXX     psMutexInit(ro->variance);
    3131    // XXX }
    3232
    … …  
    4343    // XXX     psMutexDestroy(ro->image);
    4444    // XXX }
    45     // XXX if (ro->weight) {
    46     // XXX     psMutexDestroy(ro->weight);
     45    // XXX if (ro->variance) {
     46    // XXX     psMutexDestroy(ro->variance);
    4747    // XXX }
    4848

  • trunk/psModules/src/objects/pmResiduals.c

    r21183 r21363  
    44 *
    55 * @author EAM, IfA
    6  * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2009-01-27 06:39:38 $
     6 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:25 $
    88 * Copyright 2004 IfA, University of Hawaii
    99 */
    … …  
    2323    psFree (resid->Rx);
    2424    psFree (resid->Ry);
    25     psFree (resid->weight);
     25    psFree (resid->variance);
    2626    psFree (resid->mask);
    2727    return;
    … …  
    4242    resid->Rx  = psImageAlloc (nX, nY, PS_TYPE_F32);
    4343    resid->Ry  = psImageAlloc (nX, nY, PS_TYPE_F32);
    44     resid->weight = psImageAlloc (nX, nY, PS_TYPE_F32);
     44    resid->variance = psImageAlloc (nX, nY, PS_TYPE_F32);
    4545    resid->mask   = psImageAlloc (nX, nY, PM_TYPE_RESID_MASK);
    4646
    … …  
    4949    resid->xBin = xBin;
    5050    resid->yBin = yBin;
    51     resid->xCenter = 0.5*(nX - 1); 
     51    resid->xCenter = 0.5*(nX - 1);
    5252    resid->yCenter = 0.5*(nY - 1);
    5353    return resid;

  • trunk/psModules/src/objects/pmResiduals.h

    r21183 r21363  
    44 *
    55 * @author EAM, IfA
    6  * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2009-01-27 06:39:38 $
     6 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2009-02-06 02:31:25 $
    88 * Copyright 2004 IfA, University of Hawaii
    99 */
    … …  
    1414/// @{
    1515
    16 /** residual tables for sources 
     16/** residual tables for sources
    1717 */
    1818typedef struct {
    … …  
    2020    psImage *Rx;
    2121    psImage *Ry;
    22     psImage *weight;
     22    psImage *variance;
    2323    psImage *mask;
    2424    int xBin;
    … …  
    3131bool psMemCheckResiduals(psPtr ptr);
    3232
    33 // macros to abstract the resid mask type : these values must be consistent 
     33// macros to abstract the resid mask type : these values must be consistent
    3434#define PM_TYPE_RESID_MASK PS_TYPE_U8        /**< the psElemType to use for mask image */
    3535#define PM_TYPE_RESID_MASK_DATA U8           /**< the data member to use for mask image */

  • trunk/psModules/src/objects/pmSource.c

    r21352 r21363  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:09:50 $
     8 *  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    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
    … …  
    924924        if (mode & PM_MODEL_OP_NOISE) {
    925925            // XXX need to scale by the gain...
    926             target = source->weight->data.F32;
     926            target = source->variance->data.F32;
    927927        }
    928928
    … …  
    946946    psImage *target = source->pixels;
    947947    if (mode & PM_MODEL_OP_NOISE) {
    948         target = source->weight;
     948        target = source->variance;
    949949    }
    950950

  • trunk/psModules/src/objects/pmSource.h

    r21183 r21363  
    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.28 $ $Name: not supported by cvs2svn $
     6 * @date $Date: 2009-02-06 02:31:25 $
    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

  • trunk/psModules/src/objects/pmSourceExtendedPars.c

    r20937 r21363  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2008-12-08 02:51:14 $
     8 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    7171    psFree(profile->radius);
    7272    psFree(profile->flux);
    73     psFree(profile->weight);
     73    psFree(profile->variance);
    7474    return;
    7575}
    … …  
    8282    profile->radius = NULL;
    8383    profile->flux = NULL;
    84     profile->weight = NULL;
     84    profile->variance = NULL;
    8585
    8686    return profile;

  • trunk/psModules/src/objects/pmSourceExtendedPars.h

    r15980 r21363  
    33 * @author EAM, IfA
    44 *
    5  * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    6  * @date $Date: 2008-01-02 20:39:04 $
     5 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     6 * @date $Date: 2009-02-06 02:31:25 $
    77 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    88 */
    … …  
    1717  psVector *radius;
    1818  psVector *flux;
    19   psVector *weight;
     19  psVector *variance;
    2020} pmSourceRadialProfile;
    2121

  • trunk/psModules/src/objects/pmSourceFitModel.c

    r21183 r21363  
    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.30 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    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;

  • trunk/psModules/src/objects/pmSourceFitSet.c

    r21183 r21363  
    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.15 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    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    }

  • trunk/psModules/src/objects/pmSourceMoments.c

    r21183 r21363  
    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.8 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    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;

  • trunk/psModules/src/objects/pmSourcePhotometry.c

    r21183 r21363  
    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.49 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2009-02-06 02:31:25 $
    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);

  • trunk/psModules/src/objects/pmSourceSky.c

    r21183 r21363  
    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.20 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:31:25 $
    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.