IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21484


Ignore:
Timestamp:
Feb 15, 2009, 9:51:40 AM (17 years ago)
Author:
eugene
Message:

fixes to detrend to update variance (from HEAD)

Location:
branches/eam_branch_20090208/psModules/src/detrend
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20090208/psModules/src/detrend/pmDetrendThreads.c

    r18960 r21484  
    4747
    4848    {
    49         psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_SHUTTER", 7);
     49        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_SHUTTER", 8);
    5050        task->function = &pmShutterCorrectionApplyScan_Threaded;
    5151        psThreadTaskAdd(task);
    … …  
    5454
    5555    {
    56         psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_FLAT", 9);
     56        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_FLAT", 10);
    5757        task->function = &pmFlatFieldScan_Threaded;
    5858        psThreadTaskAdd(task);

  • branches/eam_branch_20090208/psModules/src/detrend/pmFlatField.c

    r21183 r21484  
    1717    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
    1818
    19     psImage *inImage   = job->args->data[0]; // Input image
    20     psImage *inMask    = job->args->data[1]; // Input mask
    21     const psImage *flatImage = job->args->data[2]; // Flat-field image
    22     const psImage *flatMask  = job->args->data[3]; // Flat-field mask
     19    psImage *inImage = job->args->data[0]; // Input image
     20    psImage *inMask  = job->args->data[1]; // Input mask
     21    psImage *inVar   = job->args->data[2]; // Input variance
     22    const psImage *flatImage = job->args->data[3]; // Flat-field image
     23    const psImage *flatMask  = job->args->data[4]; // Flat-field mask
    2324
    24     psImageMaskType badFlat = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    25     int xOffset        = PS_SCALAR_VALUE(job->args->data[5],S32);
    26     int yOffset        = PS_SCALAR_VALUE(job->args->data[6],S32);
    27     int rowStart       = PS_SCALAR_VALUE(job->args->data[7],S32);
    28     int rowStop        = PS_SCALAR_VALUE(job->args->data[8],S32);
    29     return pmFlatFieldScan(inImage, inMask, flatImage, flatMask, badFlat,
     25    psImageMaskType badFlat = PS_SCALAR_VALUE(job->args->data[5], PS_TYPE_IMAGE_MASK_DATA);
     26    int xOffset        = PS_SCALAR_VALUE(job->args->data[6], S32);
     27    int yOffset        = PS_SCALAR_VALUE(job->args->data[7], S32);
     28    int rowStart       = PS_SCALAR_VALUE(job->args->data[8], S32);
     29    int rowStop        = PS_SCALAR_VALUE(job->args->data[9], S32);
     30
     31    return pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat,
    3032                           xOffset, yOffset, rowStart, rowStop);
    3133}
    3234
    33 // Macro for all PS types
    34 #define FLAT_DIVISION_CASE(TYPE, SPECIAL)       \
    35     case PS_TYPE_##TYPE:                        \
    36     for (int j = rowStart; j < rowStop; j++) { \
    37         for (int i = 0; i < inImage->numCols; i++) { \
    38             ps##TYPE flatValue = flatImage->data.TYPE[j + yOffset][i + xOffset]; \
    39             if (!isfinite(flatValue) || flatValue <= 0.0 || \
    40                 (flatMask && flatMask->data.PS_TYPE_IMAGE_MASK_DATA[j + yOffset][i + xOffset])) { \
    41                 if (inMask) { \
    42                     inMask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= badFlat; \
    43                 } \
    44                 inImage->data.TYPE[j][i] = SPECIAL; \
    45             } else { \
    46                 inImage->data.TYPE[j][i] /= flatValue; \
    47             } \
    48         } \
    49     } \
    50     break;
    5135
    52 bool pmFlatFieldScan(psImage *inImage, psImage *inMask, const psImage *flatImage, const psImage *flatMask,
    53                      psImageMaskType badFlat, int xOffset, int yOffset, int rowStart, int rowStop)
     36bool pmFlatFieldScan(psImage *inImage, psImage *inMask, psImage *inVar, const psImage *flatImage,
     37                     const psImage *flatMask, psImageMaskType badFlat,
     38                     int xOffset, int yOffset, int rowStart, int rowStop)
    5439{
    55     switch (inImage->type.type) {
    56         FLAT_DIVISION_CASE(U8,  0);
    57         FLAT_DIVISION_CASE(U16, 0);
    58         FLAT_DIVISION_CASE(U32, 0);
    59         FLAT_DIVISION_CASE(U64, 0);
    60         FLAT_DIVISION_CASE(S8,  0);
    61         FLAT_DIVISION_CASE(S16, 0);
    62         FLAT_DIVISION_CASE(S32, 0);
    63         FLAT_DIVISION_CASE(S64, 0);
    64         FLAT_DIVISION_CASE(F32, NAN);
    65         FLAT_DIVISION_CASE(F64, NAN);
    66     default:
    67         psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unsupported type for input image: %x\n",
    68                 inImage->type.type);
    69         return false;
     40    // Neglecting asserts, because inputs should have been checked already
     41
     42    int numCols = inImage->numCols;     // Number of columns
     43
     44    // Using i,j for image; x,y for flat.
     45    for (int j = rowStart, y = yOffset; j < rowStop; j++, y++) {
     46        for (int i = 0, x = xOffset; i < numCols; i++, x++) {
     47            float flatValue = 1.0 / flatImage->data.F32[y][x];
     48            if (!isfinite(flatValue) || flatValue <= 0.0 ||
     49                (flatMask && flatMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x])) {
     50                if (inMask) {
     51                    inMask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= badFlat;
     52                }
     53                inImage->data.F32[j][i] = NAN;
     54                if (inVar) {
     55                    inVar->data.F32[j][i] = NAN;
     56                }
     57            } else {
     58                inImage->data.F32[j][i] *= flatValue;
     59                if (inVar) {
     60                    inVar->data.F32[j][i] *= PS_SQR(flatValue);
     61                }
     62            }
     63        }
    7064    }
     65
    7166    return true;
    7267}
    … …  
    7469bool pmFlatField(pmReadout *in, const pmReadout *flat, psImageMaskType badFlat)
    7570{
    76     PS_ASSERT_PTR_NON_NULL(in, false);
    77     PS_ASSERT_PTR_NON_NULL(in->image, false);
    78     PS_ASSERT_IMAGE_NON_EMPTY(in->image, false);
    79     PS_ASSERT_PTR_NON_NULL(flat, false);
    80     PS_ASSERT_PTR_NON_NULL(flat->image, false);
    81     PS_ASSERT_IMAGE_NON_EMPTY(flat->image, false);
     71    PM_ASSERT_READOUT_NON_NULL(in, false);
     72    PM_ASSERT_READOUT_IMAGE(in, false);
     73    PM_ASSERT_READOUT_NON_NULL(flat, false);
     74    PM_ASSERT_READOUT_IMAGE(flat, false);
    8275    if (in->mask) {
    83         PS_ASSERT_IMAGE_TYPE(in->mask, PS_TYPE_IMAGE_MASK, false);
    84         PS_ASSERT_IMAGES_SIZE_EQUAL(in->mask, in->image, false);
     76        PM_ASSERT_READOUT_MASK(in, false);
    8577    }
    86     PS_ASSERT_IMAGE_TYPE(flat->image, in->image->type.type, false);
     78    if (in->variance) {
     79        PM_ASSERT_READOUT_VARIANCE(in, false);
     80    }
    8781    if (flat->mask) {
    88         PS_ASSERT_IMAGE_TYPE(flat->mask, PS_TYPE_IMAGE_MASK, false);
    89         PS_ASSERT_IMAGES_SIZE_EQUAL(flat->mask, flat->image, false);
     82        PM_ASSERT_READOUT_MASK(flat, false);
    9083    }
    9184
    9285    psImage *inImage   = in->image;     // Input image
    9386    psImage *inMask    = in->mask;      // Mask for input image
     87    psImage *inVar     = in->variance;  // Variance for input image
    9488    psImage *flatImage = flat->image;   // Flat-field image
    9589    psImage *flatMask  = flat->mask;    // Mask for flat-field image
    … …  
    146140          psArrayAdd(job->args, 1, inImage);
    147141          psArrayAdd(job->args, 1, inMask);
     142          psArrayAdd(job->args, 1, inVar);
    148143          psArrayAdd(job->args, 1, flatImage);
    149144          psArrayAdd(job->args, 1, flatMask);
    … …  
    159154          }
    160155          psFree(job);
    161       } else if (!pmFlatFieldScan(inImage, inMask, flatImage, flatMask, badFlat,
     156      } else if (!pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat,
    162157                                  xOffset, yOffset, rowStart, rowStop)) {
    163158          psError(PS_ERR_UNKNOWN, false, "Unable to flat-field image.");

  • branches/eam_branch_20090208/psModules/src/detrend/pmFlatField.h

    r21183 r21484  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2009-01-27 06:39:38 $
     7 * @version $Revision: 1.15.6.1 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-02-15 19:51:40 $
    99 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    3636    psImage *inImage,                   ///< Input image to correct
    3737    psImage *inMask,                    ///< Input mask image
     38    psImage *inVar,                     ///< Input variance image
    3839    const psImage *flatImage,           ///< Flat-field image
    3940    const psImage *flatMask,            ///< Flat-field mask
    40     psImageMaskType badFlag,                 ///< Mask value to give bad pixels
     41    psImageMaskType badFlag,            ///< Mask value to give bad pixels
    4142    int xOffset, int yOffset,           ///< Offset between input and flat-field
    4243    int rowStart, int rowStop           ///< Scan range

  • branches/eam_branch_20090208/psModules/src/detrend/pmShutterCorrection.c

    r21363 r21484  
    659659    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
    660660
    661     psImage *image        = job->args->data[0];
    662     const psImage *shutterImage = job->args->data[1];
    663     psImage *mask         = job->args->data[2];
    664 
    665     float exptime    = PS_SCALAR_VALUE(job->args->data[3],F32);
    666     psImageMaskType blank = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    667     int rowStart     = PS_SCALAR_VALUE(job->args->data[5],S32);
    668     int rowStop      = PS_SCALAR_VALUE(job->args->data[6],S32);
    669     return pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop);
    670 }
    671 
    672 bool pmShutterCorrectionApplyScan(psImage *image, const psImage *shutterImage, psImage *mask, float exptime,
     661    psImage *image = job->args->data[0];
     662    psImage *mask  = job->args->data[1];
     663    psImage *var   = job->args->data[2];
     664    const psImage *shutterImage = job->args->data[3];
     665
     666    float exptime    = PS_SCALAR_VALUE(job->args->data[4], F32);
     667    psImageMaskType blank = PS_SCALAR_VALUE(job->args->data[5], PS_TYPE_IMAGE_MASK_DATA);
     668    int rowStart     = PS_SCALAR_VALUE(job->args->data[6], S32);
     669    int rowStop      = PS_SCALAR_VALUE(job->args->data[7], S32);
     670    return pmShutterCorrectionApplyScan(image, mask, var, shutterImage, exptime, blank, rowStart, rowStop);
     671}
     672
     673bool pmShutterCorrectionApplyScan(psImage *image, psImage *mask, psImage *var,
     674                                  const psImage *shutterImage, float exptime,
    673675                                  psImageMaskType blank, int rowStart, int rowStop)
    674676{
     677    // Neglecting asserts because inputs should have been checked already
     678
     679    int numCols = image->numCols;       // Number of columns
     680
    675681    for (int y = rowStart; y < rowStop; y++) {
    676         for (int x = 0; x < image->numCols; x++) {
     682        for (int x = 0; x < numCols; x++) {
    677683            if (mask && !isfinite(shutterImage->data.F32[y][x])) {
    678684                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= blank;
    679685                image->data.F32[y][x] = NAN;
     686                if (var) {
     687                    var->data.F32[y][x] = NAN;
     688                }
    680689                continue;
    681690            }
    682             image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
     691            float correction = exptime / (exptime + shutterImage->data.F32[y][x]); // Correction factor
     692            image->data.F32[y][x] *= correction;
     693            if (var) {
     694                var->data.F32[y][x] *= PS_SQR(correction);
     695            }
    683696        }
    684697    }
    … …  
    688701bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter, psImageMaskType blank)
    689702{
    690     PS_ASSERT_PTR_NON_NULL(readout, false);
    691     PS_ASSERT_PTR_NON_NULL(shutter, false);
    692     PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
    693     PS_ASSERT_IMAGE_NON_NULL(shutter->image, false);
    694     PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
    695     PS_ASSERT_IMAGE_TYPE(shutter->image, PS_TYPE_F32, false);
     703    PM_ASSERT_READOUT_NON_NULL(readout, false);
     704    PM_ASSERT_READOUT_NON_NULL(shutter, false);
     705    PM_ASSERT_READOUT_IMAGE(readout, false);
     706    PM_ASSERT_READOUT_IMAGE(shutter, false);
    696707
    697708    psRegion region = psRegionSet(readout->col0, readout->col0 + readout->image->numCols,
    … …  
    731742    psImage *image = readout->image;    // Image to correct
    732743    psImage *mask = readout->mask;      // Corresponding mask
     744    psImage *var = readout->variance;   // Corresponding variance map
    733745
    734746    bool threaded = true;
    … …  
    766778                psThreadJob *job = psThreadJobAlloc("PSMODULES_DETREND_SHUTTER");
    767779                psArrayAdd(job->args, 1, image);
     780                psArrayAdd(job->args, 1, mask);
     781                psArrayAdd(job->args, 1, var);
    768782                psArrayAdd(job->args, 1, shutterImage);
    769                 psArrayAdd(job->args, 1, mask);
    770783                PS_ARRAY_ADD_SCALAR(job->args, exptime, PS_TYPE_F32);
    771784                PS_ARRAY_ADD_SCALAR(job->args, blank, PS_TYPE_IMAGE_MASK);
    … …  
    778791                }
    779792                psFree(job);
    780             } else if (!pmShutterCorrectionApplyScan(image, shutterImage, mask, exptime, blank,
     793            } else if (!pmShutterCorrectionApplyScan(image, mask, var, shutterImage, exptime, blank,
    781794                                                     rowStart, rowStop)) {
    782795                psError(PS_ERR_UNKNOWN, false, "Unable to apply shutter correction.");

  • branches/eam_branch_20090208/psModules/src/detrend/pmShutterCorrection.h

    r21183 r21484  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2009-01-27 06:39:38 $
     7 * @version $Revision: 1.22.6.1 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2009-02-15 19:51:40 $
    99 * Copyright 2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    9191    float offref,                       ///< Reference time offset
    9292    int nIter,                          ///< Number of iterations
    93     float rej                           ///< Rejection threshold (sigma)
     93    float rej                           ///< Rejection threshold (sigma)
    9494    );
    9595
    … …  
    131131bool pmShutterCorrectionApplyScan(
    132132    psImage *image,                     ///< Input image to correct
     133    psImage *mask,                      ///< Input mask image
     134    psImage *var,                       ///< Input variance image
    133135    const psImage *shutterImage,        ///< Shutter correction image
    134     psImage *mask,                      ///< Input mask image
    135136    float exptime,                      ///< Exposure time to which to correct
    136137    psImageMaskType blank,                   ///< Mask value to give blank pixels
Note: See TracChangeset for help on using the changeset viewer.