IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21365


Ignore:
Timestamp:
Feb 5, 2009, 4:47:00 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
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • trunk/ppMerge/src/ppMergeArguments.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    2525            "\tIMAGE(STR):     Image filename\n"
    2626            "\tMASK(STR):      Mask filename\n"
    27             "\tWEIGHT(STR)     Weight map filename\n"
    28             "where MASK and WEIGHT are optional.",
     27            "\tVARIANCE(STR)   Variance map filename\n"
     28            "where MASK and VARIANCE are optional.",
    2929            program);
    3030    fprintf(stderr, "\n");
     
    155155    psMetadataAddF32(arguments, PS_LIST_TAIL, "-frachigh", 0, "Fraction of low pixels to discard", NAN);
    156156    psMetadataAddS32(arguments, PS_LIST_TAIL, "-nkeep",    0, "Minimum number of pixels in stack to keep", 0);
    157     psMetadataAddBool(arguments, PS_LIST_TAIL, "-weights", 0, "Use image weights in combination?", false);
     157    psMetadataAddBool(arguments, PS_LIST_TAIL, "-variances", 0, "Use image variances in combination?", false);
    158158
    159159    // XXX EAM : not clear this should be allowed on the command line.
     
    300300    VALUE_ARG_RECIPE_FLOAT("-frachigh", "FRACHIGH", F32);
    301301    VALUE_ARG_RECIPE_INT("-nkeep",      "NKEEP",    S32, 0);
    302     VALUE_ARG_RECIPE_BOOL("-weights",   "WEIGHTS");
     302    VALUE_ARG_RECIPE_BOOL("-variances", "VARIANCES");
    303303
    304304    // XXX we do not supply this on the command line
  • trunk/ppMerge/src/ppMergeCamera.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    9696{
    9797    bool haveMasks = false;             // Do we have masks?
    98     bool haveWeights = false;           // Do we have weight maps?
     98    bool haveVariances = false;           // Do we have variance maps?
    9999
    100100    ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); ///< Type of frame
     
    121121
    122122        bool mdok;
    123         psString mask = psMetadataLookupStr(&mdok, input, "MASK"); ///< Name of mask
    124         psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); ///< Name of weight map
     123        psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
     124        psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
    125125
    126126        // Add the image file
     
    164164        }
    165165
    166         // Optionally add the weight file
    167         if (weight && strlen(weight) > 0) {
    168             haveWeights = true;
    169             psArray *weightFiles = psArrayAlloc(1); ///< Array of filenames for this FPA
    170             weightFiles->data[0] = psMemIncrRefCounter(weight);
    171             psMetadataAddArray(config->arguments, PS_LIST_TAIL, "WEIGHT.FILENAMES", PS_META_REPLACE,
    172                                "Filenames of weight files", weightFiles);
    173             psFree(weightFiles);
     166        // Optionally add the variance file
     167        if (variance && strlen(variance) > 0) {
     168            haveVariances = true;
     169            psArray *varianceFiles = psArrayAlloc(1); // Array of filenames for this FPA
     170            varianceFiles->data[0] = psMemIncrRefCounter(variance);
     171            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "VARIANCE.FILENAMES", PS_META_REPLACE,
     172                               "Filenames of variance files", varianceFiles);
     173            psFree(varianceFiles);
    174174
    175175            bool status;
    176             pmFPAfile *weightFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPMERGE.INPUT.WEIGHT",
    177                                                           "WEIGHT.FILENAMES");
     176            pmFPAfile *varianceFile = pmFPAfileBindFromArgs(&status, imageFile, config,
     177                                                            "PPMERGE.INPUT.VARIANCE", "VARIANCE.FILENAMES");
    178178            if (!status) {
    179                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from weight %d (%s)", numFiles, weight);
    180                 return false;
    181             }
    182             if (weightFile->type != PM_FPA_FILE_WEIGHT) {
    183                 psError(PS_ERR_IO, true, "PPMERGE.INPUT.WEIGHT is not of type WEIGHT");
    184                 return false;
    185             }
    186             haveWeights = true;
     179                psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)",
     180                        numFiles, variance);
     181                return false;
     182            }
     183            if (varianceFile->type != PM_FPA_FILE_VARIANCE) {
     184                psError(PS_ERR_IO, true, "PPMERGE.INPUT.VARIANCE is not of type VARIANCE");
     185                return false;
     186            }
     187            haveVariances = true;
    187188        }
    188189
     
    194195        psMetadataRemoveKey(config->arguments, "MASK.FILENAMES");
    195196    }
    196     if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) {
    197         psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES");
     197    if (psMetadataLookup(config->arguments, "VARIANCE.FILENAMES")) {
     198        psMetadataRemoveKey(config->arguments, "VARIANCE.FILENAMES");
    198199    }
    199200    if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) {
     
    203204    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", numFiles);
    204205    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.MASKS", 0, "Got input masks?", haveMasks);
    205     psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.WEIGHTS", 0,
    206                       "Got input weights?", haveWeights);
     206    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.VARIANCES", 0,
     207                      "Got input variances?", haveVariances);
    207208
    208209
  • trunk/ppMerge/src/ppMergeFiles.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    1313#include "ppMerge.h"
    1414
    15 const char *allFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT",
     15const char *allFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.VARIANCE",
    1616                           "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA",
    17                            NULL };      ///< All files
    18 const char *inputFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT",
    19                              NULL };    ///< Input files
     17                           NULL };      // All files
     18const char *inputFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.VARIANCE",
     19                             NULL };    // Input files
    2020const char *outputFiles[] = { "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA",
    2121                              NULL };   ///< Output files
     
    5151        }
    5252    }
    53     if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS")) {
    54         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", num);
    55         if (!pmReadoutReadChunkWeight(readout, file->fits, 0, rows, 0, config)) {
    56             psError(PS_ERR_UNKNOWN, false, "Unable to read readout weight.");
     53    if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.VARIANCES")) {
     54        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.VARIANCE", num);
     55        if (!pmReadoutReadChunkVariance(readout, file->fits, 0, rows, 0, config)) {
     56            psError(PS_ERR_UNKNOWN, false, "Unable to read readout variance.");
    5757            return false;
    5858        }
     
    9191        psFree(fileView);
    9292    }
    93     if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS")) {
    94         pmFPAfile *weight = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", num); // Weight file
    95         pmFPAview *fileView = pmFPAviewForLevel(weight->fileLevel, view);
    96         if (!pmFPAfileOpen(weight, fileView, config)) {
    97             psError(PS_ERR_UNKNOWN, false, "Unable to open weight file %d", num);
     93    if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.VARIANCES")) {
     94        pmFPAfile *variance = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.VARIANCE",
     95                                                    num); // Variance file
     96        pmFPAview *fileView = pmFPAviewForLevel(variance->fileLevel, view);
     97        if (!pmFPAfileOpen(variance, fileView, config)) {
     98            psError(PS_ERR_UNKNOWN, false, "Unable to open variance file %d", num);
    9899            psFree(fileView);
    99100            return false;
     
    127128    bool mdok;                          // Status of MD lookup
    128129    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
    129     bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Got weights?
     130    bool haveVariances = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.VARIANCES"); // Got variances?
    130131
    131132    const char **fileList = selectFiles(files); // Files to activate
     
    134135            continue;
    135136        }
    136         if (!haveWeights && strcmp(fileList[i], "PPMERGE.INPUT.WEIGHT") == 0) {
     137        if (!haveVariances && strcmp(fileList[i], "PPMERGE.INPUT.VARIANCE") == 0) {
    137138            continue;
    138139        }
     
    163164    bool mdok;                          // Status of MD lookup
    164165    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
    165     bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Got weights?
     166    bool haveVariances = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.VARIANCES"); // Got variances?
    166167
    167168    psList *list = psListAlloc(NULL);   // List of files
     
    171172            continue;
    172173        }
    173         if (!haveWeights && strcmp(fileList[i], "PPMERGE.INPUT.WEIGHT") == 0) {
     174        if (!haveVariances && strcmp(fileList[i], "PPMERGE.INPUT.VARIANCE") == 0) {
    174175            continue;
    175176        }
  • trunk/ppMerge/src/ppMergeLoop_Threaded.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    3636    bool mdok;                          ///< Status of MD lookup
    3737    bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks?
    38     bool haveWeights = psMetadataLookupBool(&mdok, arguments, "INPUTS.WEIGHTS"); ///< Do we have weights?
    39 
    40     psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); ///< Input images
    41     psArray *masks = NULL, *weights = NULL; ///< Input masks and weights
     38    bool haveVariances = psMetadataLookupBool(&mdok, arguments, "INPUTS.VARIANCES"); // Do we have variances?
     39
     40    psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images
     41    psArray *masks = NULL, *variances = NULL; // Input masks and variances
    4242    if (haveMasks) {
    4343        masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT);
    4444    }
    45     if (haveWeights) {
    46         weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT);
     45    if (haveVariances) {
     46        variances = ppMergeFileDataLevel(config, "PPMERGE.INPUT.VARIANCE", PM_FPA_LEVEL_READOUT);
    4747    }
    4848
     
    5050    if (!mdok) nThreads = 0;
    5151
    52     /** General combination parameters */
    53     int iter = psMetadataLookupS32(NULL, arguments, "ITER"); ///< Number of rejection iterations
    54     float rej = psMetadataLookupF32(NULL, arguments, "REJ"); ///< Rejection level
    55     float fraclow = psMetadataLookupF32(NULL, arguments, "FRACLOW"); ///< Reject fraction of low pixels
    56     float frachigh = psMetadataLookupF32(NULL, arguments, "FRACHIGH"); ///< Reject fraction of hi pixels
    57     int nKeep = psMetadataLookupS32(NULL, arguments, "NKEEP"); ///< Minimum number of values to keep
    58     psStatsOptions combineStat = psMetadataLookupS32(NULL, arguments, "COMBINE"); ///< Combination statistic
    59     bool useWeights = psMetadataLookupBool(NULL, arguments, "WEIGHTS"); ///< Use weights?
    60 
    61     /** Fringe parameters */
    62     int fringeNum = psMetadataLookupS32(NULL, arguments, "FRINGE.NUM"); ///< Number of fringe points
    63     int fringeSize = psMetadataLookupS32(NULL, arguments, "FRINGE.SIZE"); ///< Size of fringe regions
    64     int fringeSmoothX = psMetadataLookupS32(NULL, arguments, "FRINGE.XSMOOTH"); ///< Smoothing regions in x
    65     int fringeSmoothY = psMetadataLookupS32(NULL, arguments, "FRINGE.YSMOOTH"); ///< Smoothing regions in y
     52    // General combination parameters
     53    int iter = psMetadataLookupS32(NULL, arguments, "ITER"); // Number of rejection iterations
     54    float rej = psMetadataLookupF32(NULL, arguments, "REJ"); // Rejection level
     55    float fraclow = psMetadataLookupF32(NULL, arguments, "FRACLOW"); // Reject fraction of low pixels
     56    float frachigh = psMetadataLookupF32(NULL, arguments, "FRACHIGH"); // Reject fraction of hi pixels
     57    int nKeep = psMetadataLookupS32(NULL, arguments, "NKEEP"); // Minimum number of values to keep
     58    psStatsOptions combineStat = psMetadataLookupS32(NULL, arguments, "COMBINE"); // Combination statistic
     59    bool useVariances = psMetadataLookupBool(NULL, arguments, "VARIANCES"); // Use variances?
     60
     61    // Fringe parameters
     62    int fringeNum = psMetadataLookupS32(NULL, arguments, "FRINGE.NUM"); // Number of fringe points
     63    int fringeSize = psMetadataLookupS32(NULL, arguments, "FRINGE.SIZE"); // Size of fringe regions
     64    int fringeSmoothX = psMetadataLookupS32(NULL, arguments, "FRINGE.XSMOOTH"); // Smoothing regions in x
     65    int fringeSmoothY = psMetadataLookupS32(NULL, arguments, "FRINGE.YSMOOTH"); // Smoothing regions in y
    6666
    6767    // set the mask and mark bit values based on the named masks
    68     psImageMaskType maskVal;
    69     psImageMaskType markVal;
    70     if (!pmConfigMaskSetBits (&maskVal, &markVal, config)) {
     68    psImageMaskType maskVal, markVal;
     69    if (!pmConfigMaskSetBits(&maskVal, &markVal, config)) {
    7170        psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values");
    7271        return false;
     
    8180    combination->iter = iter;
    8281    combination->rej = rej;
    83     combination->weights = useWeights;
     82    combination->variances = useVariances;
    8483
    8584    psMetadata *stats = NULL;           ///< Statistics for output
     
    344343                psFree(countsRO);
    345344
    346                 // XXX EAM 2009.01.18 : sigmaCell and countsCell need to have their concepts copied from outCell.
    347                 // This was causing segfaults for VYSOS5; Why did this ever work for SIMTEST?
    348                 if (!pmConceptsCopyCell(countsCell, outCell)) {
    349                     psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts.");
    350                     psFree(outRO);
    351                     goto ERROR;
    352                 }
     345                // XXX EAM 2009.01.18 : sigmaCell and countsCell need to have their concepts copied from outCell.
     346                // This was causing segfaults for VYSOS5; Why did this ever work for SIMTEST?
     347                if (!pmConceptsCopyCell(countsCell, outCell)) {
     348                    psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts.");
     349                    psFree(outRO);
     350                    goto ERROR;
     351                }
    353352
    354353                pmCell *sigmaCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.SIGMA");
     
    360359                psFree(sigmaRO);
    361360
    362                 if (!pmConceptsCopyCell(sigmaCell, outCell)) {
    363                     psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts.");
    364                     psFree(outRO);
    365                     goto ERROR;
    366                 }
     361                if (!pmConceptsCopyCell(sigmaCell, outCell)) {
     362                    psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts.");
     363                    psFree(outRO);
     364                    goto ERROR;
     365                }
    367366            }
    368367
     
    447446    psFree(inputs);
    448447    psFree(masks);
    449     psFree(weights);
     448    psFree(variances);
    450449    psFree(stats);
    451450    return true;
     
    456455    psFree(inputs);
    457456    psFree(masks);
    458     psFree(weights);
     457    psFree(variances);
    459458    psFree(stats);
    460459    return false;
  • trunk/ppMerge/src/ppMergeMask.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    3838    float smoothScale = psMetadataLookupF32(&mdok, config->arguments, "MASK.SMOOTH.SCALE"); ///< Radius to grow mask
    3939
    40     psImageMaskType markVal;
    41     psImageMaskType maskValRaw;
    42     if (!pmConfigMaskSetBits (&maskValRaw, &markVal, config)) {
    43         psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values");
    44         return false;
     40    psImageMaskType markVal, maskValRaw;
     41    if (!pmConfigMaskSetBits(&maskValRaw, &markVal, config)) {
     42        psError(PS_ERR_UNKNOWN, true, "Unable to define the mask bit values");
     43        return false;
    4544    }
    4645
     
    4847    psImageMaskType maskValOut = pmConfigMaskGet (maskOutName, config);
    4948    if (!maskValOut) {
    50         psError (PS_ERR_UNKNOWN, true, "Undefined output mask bit value");
    51         return false;
    52     }
    53    
    54     psStats *statistics = psStatsAlloc(meanStat | stdevStat); ///< Statistics for background
     49        psError (PS_ERR_UNKNOWN, true, "Undefined output mask bit value");
     50        return false;
     51    }
     52
     53    psStats *statistics = psStatsAlloc(meanStat | stdevStat); // Statistics for background
    5554
    5655    psString outName = ppMergeOutputFile(config); ///< Name of output file
    5756    pmChip *outChip = pmFPAfileThisChip(config->files, view, outName); ///< Output chip
    5857    psFree(outName);
    59    
     58
    6059    int numCells = 1;
    6160    if (chipStats) {
    62         // count the number of active cells for this chip:
    63         numCells = 0;
    64         for (int i = 0; i < outChip->cells->n; i++) {
    65             pmCell *cell = outChip->cells->data[i];
    66             if (!cell->process) continue;
    67             numCells ++;
    68         }
     61        // count the number of active cells for this chip:
     62        numCells = 0;
     63        for (int i = 0; i < outChip->cells->n; i++) {
     64            pmCell *cell = outChip->cells->data[i];
     65            if (!cell->process) continue;
     66            numCells ++;
     67        }
    6968    }
    7069
     
    8584        while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
    8685
    87             // the output FPA structure carries the information about which cells to process
    88             pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); ///< Output cell
    89             if (!outCell->process) continue;
     86            // the output FPA structure carries the information about which cells to process
     87            pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell
     88            if (!outCell->process) continue;
    9089
    9190            pmHDU *hdu = pmHDUFromCell(inCell); ///< HDU for cell
     
    150149                int y = pixel / numCols;
    151150                if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValRaw)) continue;
    152                 if (outMask && (outMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValOut)) continue;
    153                 if (!isfinite(image->data.F32[y][x])) continue;
     151                if (outMask && (outMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValOut)) continue;
     152                if (!isfinite(image->data.F32[y][x])) continue;
    154153
    155154                values->data.F32[valueIndex++] = image->data.F32[y][x];
     
    165164                }
    166165
    167                 float mean = psStatsGetValue(statistics, meanStat);
    168                 float stdev = psStatsGetValue(statistics, stdevStat);
    169 
    170                 // this function increments the count for each suspect pixel in each input plane
    171                 // maskValRaw is used to test for valid input pixels
     166                float mean = psStatsGetValue(statistics, meanStat);
     167                float stdev = psStatsGetValue(statistics, stdevStat);
     168
     169                // this function increments the count for each suspect pixel in each input plane
     170                // maskValRaw is used to test for valid input pixels
    172171                if (!pmMaskFlagSuspectPixels(outRO, readout, mean, stdev, maskSuspect, maskValRaw)) {
    173172                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
     
    197196            inView->cell = -1;
    198197            while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
    199 
    200                 // the output FPA structure carries the information about which cells to process
    201                 pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); ///< Output cell
    202                 if (!outCell->process) continue;
     198                // the output FPA structure carries the information about which cells to process
     199                pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell
     200                if (!outCell->process) continue;
    203201
    204202                pmHDU *hdu = pmHDUFromCell(inCell); ///< HDU for cell
     
    210208                pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK");
    211209
    212                 float mean = psStatsGetValue(statistics, meanStat);
    213                 float stdev = psStatsGetValue(statistics, stdevStat);
    214 
    215                 if (!pmMaskFlagSuspectPixels(outRO, readout, mean, stdev, maskSuspect, maskValRaw)) {
     210                float mean = psStatsGetValue(statistics, meanStat);
     211                float stdev = psStatsGetValue(statistics, stdevStat);
     212
     213                if (!pmMaskFlagSuspectPixels(outRO, readout, mean, stdev, maskSuspect, maskValRaw)) {
    216214                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
    217215                    goto MERGE_MASK_ERROR;
     
    243241    while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) {
    244242
    245         // skip inactive cells
    246         if (!outCell->process) continue;
     243        // skip inactive cells
     244        if (!outCell->process) continue;
    247245
    248246        pmHDU *hdu = pmHDUFromCell(outCell); ///< HDU for cell
     
    254252        pmReadout *outRO = outCell->readouts->data[0]; ///< Output readout
    255253
    256         if (smoothSuspect) {
    257             // XXX test output of suspect pixel image
    258             psImage *suspects = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); ///< Suspect img
    259             assert (suspects);
    260             psImageSmooth (suspects, smoothScale, 3); ///< extend smoothing region to 3-sigma
    261         }
    262 
    263         // set the bad pixels to the value 'maskVal'
     254        if (smoothSuspect) {
     255            // XXX test output of suspect pixel image
     256            psImage *suspects = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img
     257            assert (suspects);
     258            psImageSmooth (suspects, smoothScale, 3); // extend smoothing region to 3-sigma
     259        }
     260
     261        // set the bad pixels to the value 'maskVal'
    264262        if (!pmMaskIdentifyBadPixels(outRO, maskValOut, maskBad, maskMode)) {
    265263            psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels");
     
    352350    assert(config);
    353351
    354     bool mdok;                          ///< Status of MD lookup
    355     int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); ///< Number of inputs
    356     bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); ///< Do we have masks?
    357     bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); ///< Do we have weights?
    358     int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); ///< Number of rejection iterations
     352    bool mdok;                          // Status of MD lookup
     353    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     354    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
     355    bool haveVariances = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.VARIANCES"); // Got variances?
     356    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations
    359357
    360358    PS_ASSERT_INT_POSITIVE(iter, false);
     
    375373        psFree(masks);
    376374    }
    377     if (haveWeights) {
    378         psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT);
    379         psFree(weights);
     375    if (haveVariances) {
     376        psArray *variances = ppMergeFileDataLevel(config, "PPMERGE.INPUT.VARIANCE", PM_FPA_LEVEL_READOUT);
     377        psFree(variances);
    380378    }
    381379
     
    417415        }
    418416
    419         if (outChip->data_exists) {
    420             psList *inChips = psListAlloc(NULL);
    421             for (int i=0; i < numFiles; i++) {
    422                 pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); ///< Input file
    423                 pmChip *chip = pmFPAviewThisChip(view, file->fpa);
    424                 psListAdd(inChips, PS_LIST_TAIL, chip);
    425             }
    426 
    427             // XXX I need to call pmConfigMaskWriteHeader for the PHU somewhere, after it is created!
    428 
    429             if (!pmConceptsAverageChips(outChip, inChips, true)) {
    430                 psError(PS_ERR_UNKNOWN, false, "Unable to average Chip concepts.");
    431                 psFree(inChips);
    432                 goto PPMERGE_MASK_ERROR;
    433             }
    434             psFree(inChips);
    435         }
     417        if (outChip->data_exists) {
     418            psList *inChips = psListAlloc(NULL);
     419            for (int i=0; i < numFiles; i++) {
     420                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
     421                pmChip *chip = pmFPAviewThisChip(view, file->fpa);
     422                psListAdd(inChips, PS_LIST_TAIL, chip);
     423            }
     424
     425            // XXX I need to call pmConfigMaskWriteHeader for the PHU somewhere, after it is created!
     426
     427            if (!pmConceptsAverageChips(outChip, inChips, true)) {
     428                psError(PS_ERR_UNKNOWN, false, "Unable to average Chip concepts.");
     429                psFree(inChips);
     430                goto PPMERGE_MASK_ERROR;
     431            }
     432            psFree(inChips);
     433        }
    436434        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    437435            goto PPMERGE_MASK_ERROR;
  • trunk/ppMerge/src/ppMergeReadChunk.c

    r21244 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-01 21:43:05 $
     8 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    2020
    2121    bool mdok;
    22     bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); ///< Do we have masks?
    23     bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS");///< Do we have weights?
    24     int rows = psMetadataLookupS32(NULL, config->arguments, "ROWS"); ///< Number of rows to read per chunk
     22    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
     23    bool haveVariances = psMetadataLookupBool(&mdok, config->arguments,
     24                                              "INPUTS.VARIANCES");// Do we have variances?
     25    int rows = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of rows to read per chunk
    2526
    2627    // select an available fileGroup
     
    5253                // override the recorded last scan
    5354                inRO->thisImageScan  = fileGroup->firstScan;
    54                 inRO->thisWeightScan = fileGroup->firstScan;
     55                inRO->thisVarianceScan = fileGroup->firstScan;
    5556                inRO->thisMaskScan   = fileGroup->firstScan;
    5657                inRO->forceScan      = true;
     
    8182                }
    8283
    83                 if (haveWeights && pmReadoutMoreWeight(inRO, file->fits, 0, rows, config)) {
     84                if (haveVariances && pmReadoutMoreVariance(inRO, file->fits, 0, rows, config)) {
    8485                    keepReading = true;
    85                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i);
    86                     if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) {
    87                         psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d",
     86                    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.VARIANCE", i);
     87                    if (!pmReadoutReadChunkVariance(inRO, file->fits, 0, rows, 0, config)) {
     88                        psError(PS_ERR_IO, false,
     89                                "Unable to read chunk %d for file PPMERGE.INPUT.VARIANCE %d",
    8890                                numChunk, i);
    8991                        *status = false;
  • trunk/ppMerge/src/ppMergeScaleZero.c

    r21265 r21365  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-02 20:26:26 $
     8 *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-06 02:44:31 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     
    5656    for (int i = 0; i < numInputs; i++) {
    5757        pmFPAfileActivate(config->files, false, NULL);
    58         psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); ///< Activated files
    59         pmFPAfile *input = files->data[0]; ///< Representative file; should be the image (not mask or weight)
    60         pmFPA *fpa = input->fpa;        ///< FPA of interest
    61         view = pmFPAviewAlloc(0);       ///< View to component of interest
    62         int cellNum = 0;                ///< Index for cell
     58        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Activated files
     59        pmFPAfile *input = files->data[0]; // Representative file; should be the image (not mask or variance)
     60        pmFPA *fpa = input->fpa;        // FPA of interest
     61        view = pmFPAviewAlloc(0);       // View to component of interest
     62        int cellNum = 0;                // Index for cell
    6363        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    6464            goto ERROR;
  • trunk/ppSim/src/ppSimInsertStars.c

    r20367 r21365  
    3939    float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate
    4040    if (isnan(skyRate)) {
    41         float zp      = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok);
    42         float scale   = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE");     assert (mdok);
    43         float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
     41        float zp      = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok);
     42        float scale   = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE");     assert (mdok);
     43        float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
    4444        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    4545    }
     
    7575    FILE *outfile = fopen (outname, "w");
    7676
    77     // add sources to the readout image & weight
     77    // add sources to the readout image & variance
    7878    psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
    7979    for (long i = 0; i < stars->n; i++) {
     
    118118        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    119119
    120         // this value is the pure (input) flux, and is saved in the output source cmf files
     120        // this value is the pure (input) flux, and is saved in the output source cmf files
    121121        source->psfMag = -2.5*log10(star->flux);
    122122        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
     
    132132        // Blow away the image parts of the source, which makes the memory explode
    133133        RESET(source->pixels);
    134         RESET(source->weight);
     134        RESET(source->variance);
    135135        RESET(source->maskObj);
    136136        RESET(source->maskView);
  • trunk/ppSim/src/ppSimLoadForceSources.c

    r18011 r21365  
    5656    // Select the spots within range of this readout.  Project the spots to this chip
    5757    for (int i = 0; i < spots->n; i++) {
    58         pmAstromObj *obj = spots->data[i];
     58        pmAstromObj *obj = spots->data[i];
    5959
    60         psProject (obj->TP, obj->sky, fpa->toSky);
    61         psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP);
    62         psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP);
     60        psProject (obj->TP, obj->sky, fpa->toSky);
     61        psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP);
     62        psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP);
    6363
    64         // limit the X,Y range of the objs to the selected chip
    65         if (obj->chip->x < minX) continue;
    66         if (obj->chip->x > maxX) continue;
    67         if (obj->chip->y < minY) continue;
    68         if (obj->chip->y > maxY) continue;
     64        // limit the X,Y range of the objs to the selected chip
     65        if (obj->chip->x < minX) continue;
     66        if (obj->chip->x > maxX) continue;
     67        if (obj->chip->y < minY) continue;
     68        if (obj->chip->y > maxY) continue;
    6969
    70         // convert the pmAstromObj to pmSource
     70        // convert the pmAstromObj to pmSource
    7171
    7272        // instantiate a model for the PSF at this location, Io = 1.0
     
    7474
    7575        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
    76         // XXX use a fixed radius??
     76        // XXX use a fixed radius??
    7777        // float radius = model->modelRadius (model->params, roughNoise);
    7878        // radius = PS_MAX (radius, 1.0);
    79         float radius = 5.0;
     79        float radius = 5.0;
    8080
    8181        // construct a source, with model flux pixels set, based on the model
     
    8787        // psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    8888
    89         // these are not really needed since we will be fitting for them
     89        // these are not really needed since we will be fitting for them
    9090        source->psfMag = NAN;
    9191        source->errMag = NAN;
    9292
    9393        // insert the source flux in the image
    94         // XXX not sure the offset is really 0,0
     94        // XXX not sure the offset is really 0,0
    9595        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0);
    9696        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0);
     
    9898        // Blow away the image parts of the source, which makes the memory explode
    9999        RESET(source->pixels);
    100         RESET(source->weight);
     100        RESET(source->variance);
    101101        RESET(source->maskObj);
    102102        RESET(source->maskView);
     
    105105        RESET(source->blends);
    106106
    107         psArrayAdd (sources, 100, source);
     107        psArrayAdd (sources, 100, source);
    108108        psFree(source);                 // Drop local reference
    109109    }
  • trunk/ppSim/src/ppSimLoop.c

    r18011 r21365  
    3434    psArray *galaxies = psArrayAllocEmpty (1);
    3535    if (type == PPSIM_TYPE_OBJECT) {
    36         // Load forced-photometry positions (these are placed on fpa->analysis for use in ppSimPhotomReadout)
    37         if (!ppSimLoadSpots (fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots");
    38 
    39         // Load catalogue stars
     36        // Load forced-photometry positions (these are placed on fpa->analysis for use in ppSimPhotomReadout)
     37        if (!ppSimLoadSpots (fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots");
     38
     39        // Load catalogue stars
    4040        if (!ppSimLoadStars (stars, fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load catalog stars");
    4141
    42         // Add random stars
     42        // Add random stars
    4343        if (!ppSimMakeStars (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random stars");
    4444
    45         // Add random galaxies
     45        // Add random galaxies
    4646        if (!ppSimMakeGalaxies (galaxies, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random galaxies");
    4747    }
     
    105105                // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately
    106106                readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels
    107                 readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
     107                readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
    108108
    109109                psImageInit (readout->image, 0.0);
    110                 psImageInit (readout->weight, 0.0);
     110                psImageInit (readout->variance, 0.0);
    111111
    112112                psFree(readout);        // Drop reference
     
    115115            psVector *biasCols = ppSimMakeBiassec (cell, config);
    116116
    117             pmReadout *readout;
    118             while ((readout = pmFPAviewNextReadout (view, fpa, 1))) {
    119 
    120                 // if we have not read in a weight or generated a fake image above, we need to
     117            pmReadout *readout;
     118            while ((readout = pmFPAviewNextReadout (view, fpa, 1))) {
     119
     120                // if we have not read in a variance or generated a fake image above, we need to
    121121                // build one here
    122                 if (!readout->weight) {
    123                     if (!pmReadoutGenerateWeight(readout, true)) {
    124                         psError (PS_ERR_UNKNOWN, false, "trouble creating weight");
     122                if (!readout->variance) {
     123                    if (!pmReadoutGenerateVariance(readout, true)) {
     124                        psError (PS_ERR_UNKNOWN, false, "trouble creating variance");
    125125                        return false;
    126126                    }
     
    133133
    134134                psVector *biasRows = ppSimMakeBias (&status, readout, config, rng);
    135                 if (!status) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure");
     135                if (!status) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure");
    136136                if (type == PPSIM_TYPE_BIAS) goto done;
    137137
    138                 if (!ppSimMakeDark (readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure");
     138                if (!ppSimMakeDark (readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure");
    139139                if (type == PPSIM_TYPE_DARK) goto done;
    140140
     
    153153
    154154            done:
    155                 if (!ppSimAddNoise(readout->image, readout->weight, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise");
     155                if (!ppSimAddNoise(readout->image, readout->variance, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise");
    156156                if (!ppSimSaturate(readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem setting saturation levels");
    157157
     
    165165                readout->parent->parent->data_exists = true;
    166166
    167                 // if there is an input image, merge it with the simulated image
    168                 if (!ppSimMergeReadouts (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image");
     167                // if there is an input image, merge it with the simulated image
     168                if (!ppSimMergeReadouts (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image");
    169169            }
    170170            psFree(biasCols);
     
    173173
    174174            if (cell->hdu) {
    175                 // XXX only do this if there is no INPUT image?
     175                // XXX only do this if there is no INPUT image?
    176176                if (!ppSimInitHeader(config, NULL, NULL, cell)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header");
    177177            }
     
    186186        }
    187187
    188         // XXX why no UpdateConceptsChip??
     188        // XXX why no UpdateConceptsChip??
    189189
    190190        if (chip->hdu) {
    191             // XXX only do this if there is no INPUT image
     191            // XXX only do this if there is no INPUT image
    192192            if (!ppSimInitHeader(config, NULL, chip, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header");
    193193        }
    194194
    195195        // we perform photometry on the readouts of this chip in the output
    196         if (!ppSimPhotom (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem performing photometry");
     196        if (!ppSimPhotom (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem performing photometry");
    197197
    198198        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
  • trunk/ppSim/src/ppSimMakeBias.c

    r18011 r21365  
    2020    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    2121    if (isnan(readnoise)) {
    22         psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
    23         readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
    24         if (!mdok) {
    25             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    26             *status = false;
    27             return NULL;
    28         }
     22        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
     23        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
     24        if (!mdok) {
     25            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
     26            *status = false;
     27            return NULL;
     28        }
    2929    }
    3030
    3131    psImage *signal = readout->image;
    32     psImage *variance = readout->weight;
     32    psImage *variance = readout->variance;
    3333
    3434    int numRows = signal->numRows;
     
    3838    psPolynomial1D *biasPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, biasOrder);
    3939    for (int j = 0; j < biasOrder + 1; j++) {
    40         biasPoly->coeff[j] = biasRange * psRandomGaussian(rng);
     40        biasPoly->coeff[j] = biasRange * psRandomGaussian(rng);
    4141    }
    4242
     
    4545
    4646    for (int y = 0; y < numRows; y++) {
    47         // Adjust bias level for this row
    48         biasRows->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) /
    49                                                   (float)numRows - 0.5) + biasLevel;
     47        // Adjust bias level for this row
     48        biasRows->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) /
     49                                                  (float)numRows - 0.5) + biasLevel;
    5050
    51         for (int x = 0; x < numCols; x++) {
     51        for (int x = 0; x < numCols; x++) {
    5252
    53             // Bias level
    54             signal->data.F32[y][x] += biasRows->data.F32[y];
    55             variance->data.F32[y][x] += PS_SQR(readnoise);
     53            // Bias level
     54            signal->data.F32[y][x] += biasRows->data.F32[y];
     55            variance->data.F32[y][x] += PS_SQR(readnoise);
    5656
    57         }
     57        }
    5858    }
    5959    psFree(biasPoly);
  • trunk/ppSim/src/ppSimMakeDark.c

    r18011 r21365  
    77
    88    psImage *signal = readout->image;
    9     psImage *variance = readout->weight;
     9    psImage *variance = readout->variance;
    1010
    1111    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     
    2929
    3030    for (int y = 0; y < signal->numRows; y++) {
    31         for (int x = 0; x < signal->numCols; x++) {
    32                        
    33             // Dark current
    34             float darkCurrent = darkRate * expTime; // Dark current accumulated
    35             signal->data.F32[y][x] += darkCurrent;
    36             variance->data.F32[y][x] += darkCurrent;
    37         }
     31        for (int x = 0; x < signal->numCols; x++) {
     32
     33            // Dark current
     34            float darkCurrent = darkRate * expTime; // Dark current accumulated
     35            signal->data.F32[y][x] += darkCurrent;
     36            variance->data.F32[y][x] += darkCurrent;
     37        }
    3838    }
    3939    return true;
  • trunk/ppSim/src/ppSimMakeSky.c

    r20366 r21365  
    99
    1010    psImage *signal = readout->image;
    11     psImage *variance = readout->weight;
     11    psImage *variance = readout->variance;
    1212
    1313    pmCell *cell = readout->parent;
     
    1919    bool sky  = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux?
    2020    bool flat = psMetadataLookupBool(&status, recipe, "FLAT"); // Apply flat-field term?
    21  
     21
    2222    float expTime      = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time
    2323
     
    2929    float skyMags      = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
    3030    if (!isnan(skyMags)) {
    31         float zp       = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status);
    32         float scale    = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status);
    33         skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
     31        float zp       = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status);
     32        float scale    = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status);
     33        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    3434    }
    3535    if (type == PPSIM_TYPE_FLAT) {
     
    3737    }
    3838
    39     int x0Chip        = psMetadataLookupS32(&status, chip->concepts, "CHIP.X0");
    40     int y0Chip        = psMetadataLookupS32(&status, chip->concepts, "CHIP.Y0");
     39    int x0Chip        = psMetadataLookupS32(&status, chip->concepts, "CHIP.X0");
     40    int y0Chip        = psMetadataLookupS32(&status, chip->concepts, "CHIP.Y0");
    4141    int xParityChip   = psMetadataLookupS32(&status, chip->concepts, "CHIP.XPARITY");
    4242    int yParityChip   = psMetadataLookupS32(&status, chip->concepts, "CHIP.YPARITY");
    4343
    44     int x0Cell        = psMetadataLookupS32(&status, cell->concepts, "CELL.X0");
    45     int y0Cell        = psMetadataLookupS32(&status, cell->concepts, "CELL.Y0");
     44    int x0Cell        = psMetadataLookupS32(&status, cell->concepts, "CELL.X0");
     45    int y0Cell        = psMetadataLookupS32(&status, cell->concepts, "CELL.Y0");
    4646    int xParityCell   = psMetadataLookupS32(&status, cell->concepts, "CELL.XPARITY");
    4747    int yParityCell   = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY");
     
    7272            // Gaussian flat-field over the FPA with flatValue = 1.0 at the field center
    7373            float flatValue = 1.0;
    74             if (flat) {
    75                 // we make the flat-field have a response of 1.0 at the field center (like a vignetting)
    76                 flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA)));
    77             }
     74            if (flat) {
     75                // we make the flat-field have a response of 1.0 at the field center (like a vignetting)
     76                flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA)));
     77            }
    7878
    79             float scatterRate = 0.0;
     79            float scatterRate = 0.0;
    8080
    81             if (sky) {
    82               // add a scattered light term to the flat-field images
    83               if (type == PPSIM_TYPE_FLAT) {
    84                   scatterRate = scatterFrac * PS_SQR(xFPA);
    85               }
     81            if (sky) {
     82              // add a scattered light term to the flat-field images
     83              if (type == PPSIM_TYPE_FLAT) {
     84                  scatterRate = scatterFrac * PS_SQR(xFPA);
     85              }
    8686
    87               // Sky background
    88               float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky
    89               signal->data.F32[y][x] += skyFlux;
    90               variance->data.F32[y][x] += skyFlux;
    91             }
     87              // Sky background
     88              float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky
     89              signal->data.F32[y][x] += skyFlux;
     90              variance->data.F32[y][x] += skyFlux;
     91            }
    9292
    93             // used later to modify the star and galaxy photometry
    94             if (expCorr) {
    95                 // exposure correction is (effective exposure time) * (flatValue)
    96               expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
    97             }
     93            // used later to modify the star and galaxy photometry
     94            if (expCorr) {
     95                // exposure correction is (effective exposure time) * (flatValue)
     96              expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
     97            }
    9898
    9999            // TO DO: Add fringes
  • trunk/ppSim/src/ppSimMergeReadouts.c

    r18011 r21365  
    1414    if (!inReadout) return true;
    1515
    16     if (!inReadout->weight) {
    17       if (!pmReadoutGenerateWeight(inReadout, true)) {
    18         psError (PS_ERR_UNKNOWN, false, "trouble creating weight");
    19         return false;
     16    if (!inReadout->variance) {
     17      if (!pmReadoutGenerateVariance(inReadout, true)) {
     18        psError (PS_ERR_UNKNOWN, false, "trouble creating variance");
     19        return false;
    2020      }
    2121    }
     
    3030
    3131    psImage *inSignal = inReadout->image;
    32     psImage *inVariance = inReadout->weight;
     32    psImage *inVariance = inReadout->variance;
    3333
    3434    psImage *outSignal = outReadout->image;
    35     psImage *outVariance = outReadout->weight;
     35    psImage *outVariance = outReadout->variance;
    3636
    3737    assert (inSignal->numRows == outSignal->numRows);
     
    3939
    4040    for (int y = 0; y < inSignal->numRows; y++) {
    41         for (int x = 0; x < inSignal->numCols; x++) {
    42             outSignal->data.F32[y][x] += inSignal->data.F32[y][x];
    43             outVariance->data.F32[y][x] += inVariance->data.F32[y][x];
    44         }
     41        for (int x = 0; x < inSignal->numCols; x++) {
     42            outSignal->data.F32[y][x] += inSignal->data.F32[y][x];
     43            outVariance->data.F32[y][x] += inVariance->data.F32[y][x];
     44        }
    4545    }
    4646    return true;
  • trunk/ppSim/src/ppSimPhotomReadoutFake.c

    r21183 r21365  
    1212    }
    1313
    14 # if 0   
     14# if 0
    1515    // set the photcode for this image
    1616    if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) {
     
    2222    // *** in this section, perform the photometry for real + fake sources on PPSIM.FAKE.CHIP ***
    2323
    24     // find the currently selected readout. 
     24    // find the currently selected readout.
    2525    // we always perform photometry on the mosaiced chip
    2626    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.CHIP");
     
    3535    PS_ASSERT_PTR_NON_NULL (injectedSources, false);
    3636
    37     // Generate the mask and weight images, including the user-defined analysis region of interest
    38     psphotSetMaskAndWeight (config, readout, recipe);
     37    // Generate the mask and variance images, including the user-defined analysis region of interest
     38    psphotSetMaskAndVariance (config, readout, recipe);
    3939
    4040    // XXX need to define the source pixels
     
    5353    psArray *realSources = psArrayAlloc (realMeasuredSources->n);
    5454    for (int i = 0; i < realMeasuredSources->n; i++) {
    55         realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
     55        realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
    5656    }
    5757
    5858    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    59     // this function uses PSPHOT.PSF.LOAD as the pmFPAfile 
     59    // this function uses PSPHOT.PSF.LOAD as the pmFPAfile
    6060    pmPSF *psf = psphotLoadPSF (config, view, recipe);
    6161    assert (psf);
    6262
    63     // remove all sources 
     63    // remove all sources
    6464    psphotRemoveAllSources (realSources, recipe);
    6565
     
    7676
    7777    // XXX fake sources should measure peak->x,y, force sources should not
    78     psImageMaskType maskVal = 0xff; 
     78    psImageMaskType maskVal = 0xff;
    7979    psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal);
    8080    ppSimDetections (significance, recipe, fakeSources);
     
    9090    psphotGuessModels (config, readout, realSources, psf);
    9191    psphotGuessModels (config, readout, fakeSources, psf);
    92    
     92
    9393    // linear fit to real + fake sources
    9494    psArray *sources = ppSimMergeSources (realSources, fakeSources);
     
    113113    pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES");
    114114    if (!fakeReadout) {
    115         fakeReadout = pmReadoutAlloc (fakeCell);
    116         psFree (fakeReadout); // there is a copy on 'cell' as well
     115        fakeReadout = pmReadoutAlloc (fakeCell);
     116        psFree (fakeReadout); // there is a copy on 'cell' as well
    117117    }
    118118    psAssert (fakeReadout, "no fakeReadout?");
  • trunk/ppSim/src/ppSimPhotomReadoutForce.c

    r21183 r21365  
    1212    }
    1313
    14 # if 0   
     14# if 0
    1515    // set the photcode for this image
    1616    if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) {
     
    2222    // *** in this section, perform the photometry for real + force sources on PPSIM.FORCE.CHIP ***
    2323
    24     // find the currently selected readout. 
     24    // find the currently selected readout.
    2525    // we always perform photometry on the mosaiced chip
    2626    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP");
     
    3535    psAssert (forceSources, "failed to load force photometry sources");
    3636
    37     // Generate the mask and weight images, including the user-defined analysis region of interest
    38     psphotSetMaskAndWeight (config, readout, recipe);
     37    // Generate the mask and variance images, including the user-defined analysis region of interest
     38    psphotSetMaskAndVariance (config, readout, recipe);
    3939
    4040    // XXX need to define the source pixels
     
    4646    psArray *realSources = psArrayAlloc (realMeasuredSources->n);
    4747    for (int i = 0; i < realMeasuredSources->n; i++) {
    48         realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
     48        realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
    4949    }
    5050
    5151    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    52     // this function uses PSPHOT.PSF.LOAD as the pmFPAfile 
     52    // this function uses PSPHOT.PSF.LOAD as the pmFPAfile
    5353    pmPSF *psf = psphotLoadPSF (config, view, recipe);
    5454    assert (psf);
    5555
    56     // remove all sources 
     56    // remove all sources
    5757    psphotRemoveAllSources (realSources, recipe);
    5858
     
    6969
    7070    // XXX fake sources should measure peak->x,y, force sources should not
    71     psImageMaskType maskVal = 0xff; 
     71    psImageMaskType maskVal = 0xff;
    7272    psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal);
    7373    ppSimDetections (significance, recipe, forceSources);
     
    8383    psphotGuessModels (config, readout, realSources, psf);
    8484    psphotGuessModels (config, readout, forceSources, psf);
    85    
     85
    8686    // linear fit to real + force sources
    8787    psArray *sources = ppSimMergeSources (realSources, forceSources);
     
    106106    pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES");
    107107    if (!forceReadout) {
    108         forceReadout = pmReadoutAlloc (forceCell);
    109         psFree (forceReadout); // there is a copy on 'cell' as well
     108        forceReadout = pmReadoutAlloc (forceCell);
     109        psFree (forceReadout); // there is a copy on 'cell' as well
    110110    }
    111111    psAssert (forceReadout, "no forceReadout?");
Note: See TracChangeset for help on using the changeset viewer.