IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 7, 2009, 4:08:25 PM (17 years ago)
Author:
Paul Price
Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:
branches/pap
Files:
82 edited
10 copied

Legend:

Unmodified
Added
Removed

  • branches/pap

  • branches/pap/psModules

    • Property svn:mergeinfo deleted

  • branches/pap/psModules/src/astrom/pmAstrometryDistortion.c

    r23487 r25027  
    151151
    152152            // also measure the L and M median positions as a representative coordinate
    153             psVectorStats (stats, L, NULL, NULL, 0);
     153            if (!psVectorStats (stats, L, NULL, NULL, 0)) {
     154                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     155                goto skip;
     156            }
    154157            grad->FP.x = stats->sampleMedian;
    155158
    156             psVectorStats (stats, M, NULL, NULL, 0);
     159            if (!psVectorStats (stats, M, NULL, NULL, 0)) {
     160                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     161                goto skip;
     162            }
    157163            grad->FP.y = stats->sampleMedian;
    158164

  • branches/pap/psModules/src/astrom/pmAstrometryObjects.c

    r21422 r25027  
    8383{
    8484    psArray *matches = psArrayAllocEmpty(x1->n);
     85    psVector *found1 = psVectorAlloc(x1->n, PS_TYPE_S8);
     86    psVector *found2 = psVectorAlloc(x2->n, PS_TYPE_S8);
    8587
    8688    const double RADIUS_SQR = PS_SQR(RADIUS);
    8789    double dX, dY, dR;
     90
     91    psVectorInit (found1, 0);
     92    psVectorInit (found2, 0);
    8893
    8994    int jStart;
    … …  
    100105        }
    101106
     107        if (found1->data.S8[i]) {
     108            i++;
     109            continue;
     110        }
     111        if (found2->data.S8[j]) {
     112            j++;
     113            continue;
     114        }
     115
    102116        jStart = j;
    103         while (fabs(dX) < RADIUS && j < x2->n) {
     117        while ((fabs(dX) < RADIUS) && (j < x2->n)) {
    104118
    105119            dX = x1->data.F64[i] - x2->data.F64[j];
    … …  
    111125                continue;
    112126            }
     127            if (found2->data.S8[j]) {
     128                j++;
     129                continue;
     130            }
    113131
    114132            // got a match; add to output list
    … …  
    117135            psFree (match);
    118136
     137            found1->data.S8[i] = 1;
     138            found2->data.S8[j] = 1;
     139
    119140            j++;
    120141        }
    … …  
    122143        i++;
    123144    }
     145    psFree (found1);
     146    psFree (found2);
     147
    124148    return (matches);
    125149}
    … …  
    420444    obj->sky  = psSphereAlloc();
    421445    obj->Mag  = 0;
     446    obj->Color= 0;
    422447    obj->dMag = 0;
    423448
    … …  
    447472    *obj->sky  = *old->sky;
    448473    obj->Mag   =  old->Mag;
     474    obj->Color =  old->Color;
    449475    obj->dMag  =  old->dMag;
    450476

  • branches/pap/psModules/src/astrom/pmAstrometryObjects.h

    r20801 r25027  
    3333typedef struct
    3434{
    35     psPlane *pix;   ///< the position in the pmReadout frame
    36     psPlane *cell;   ///< the position in the pmCell frame
    37     psPlane *chip;   ///< the position in the pmChip frame
    38     psPlane *FP;   ///< the position in the pmFPA frame
    39     psPlane *TP;   ///< the position in the tangent plane
    40     psSphere *sky;        ///< the position on the Celestial Sphere.
    41     double Mag;                         ///< object magnitude XXX what filter?
    42     double dMag;                        ///< error on object magnitude
     35    psPlane *pix;                       ///< the position in the pmReadout frame
     36    psPlane *cell;                      ///< the position in the pmCell frame
     37    psPlane *chip;                      ///< the position in the pmChip frame
     38    psPlane *FP;                        ///< the position in the pmFPA frame
     39    psPlane *TP;                        ///< the position in the tangent plane
     40    psSphere *sky;                      ///< the position on the Celestial Sphere.
     41    float Mag;                          ///< object magnitude in extracted filter
     42    float Color;                        ///< object color
     43    float dMag;                         ///< error on object magnitude
    4344}
    4445pmAstromObj;

  • branches/pap/psModules/src/astrom/pmAstrometryUtils.c

    r19991 r25027  
    4646
    4747        /* this loop uses the Newton-Raphson method to solve for Xo,Yo
    48         * it needs the high order terms to be small
    49         * Xo,Yo are in pixels;
    50         */
     48        * it needs the high order terms to be small
     49        * Xo,Yo are in pixels;
     50        */
    5151        double dPos = tol + 1;
    5252        for (int i = 0; (dPos > tol) && (i < 20); i++) {
    … …  
    6060            Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);
    6161
    62             if (!psMatrixGJSolveF32 (Alpha, Beta)) {
    63               psError(PS_ERR_UNKNOWN, false, "Unable to solve for center.");
    64               psFree (Alpha);
    65               psFree (Beta);
    66               psFree (XdX);
    67               psFree (XdY);
    68               psFree (YdX);
    69               psFree (YdY);
     62            if (!psMatrixGJSolve (Alpha, Beta)) {
     63                psError(PS_ERR_UNKNOWN, false, "Unable to solve for center.");
     64                psFree (Alpha);
     65                psFree (Beta);
     66                psFree (XdX);
     67                psFree (XdY);
     68                psFree (YdX);
     69                psFree (YdY);
    7070                return NULL;
    7171            }

  • branches/pap/psModules/src/astrom/pmAstrometryVisual.c

    r23487 r25027  
    1717#include "pmFPAExtent.h"
    1818
    19 # if (HAVE_KAPA)
    20 # include <kapa.h>
    21 # include "pmKapaPlots.h"
     19#if (HAVE_KAPA)
     20#include <kapa.h>
     21#include "pmKapaPlots.h"
    2222#include "pmVisual.h"
    2323

  • branches/pap/psModules/src/astrom/pmAstrometryWCS.c

    r21430 r25027  
    9696
    9797    if (!status1 || !status2) {
     98        Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1");
     99        Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2");
     100    }
     101
     102    if (!status1 || !status2) {
    98103        Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
    99104        Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
    … …  
    143148    int Nx = region->x1 - region->x0;
    144149    int Ny = region->y1 - region->y0;
    145     psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Mosaic Dimensions", Nx);
    146     psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Mosaic Dimensions", Ny);
     150    psMetadataAddS32 (header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "Mosaic Dimensions", Nx);
     151    psMetadataAddS32 (header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "Mosaic Dimensions", Ny);
    147152
    148153    pmAstromWCStoHeader (header, wcs);

  • branches/pap/psModules/src/astrom/pmAstrometryWCS.h

    r21159 r25027  
    6868bool psPlaneDistortIsDiagonal (psPlaneDistort *distort);
    6969
    70 # define PM_DEG_RAD 57.295779513082322
    71 # define PM_RAD_DEG  0.017453292519943
     70// XXX probably should remove these and just use the PS_ version in the code
     71# define PM_DEG_RAD PS_DEG_RAD
     72# define PM_RAD_DEG PS_RAD_DEG
    7273
    7374/// @}

  • branches/pap/psModules/src/camera/Makefile.am

    r19695 r25027  
    2323        pmFPAfileIO.c \
    2424        pmFPAfileFitsIO.c \
     25        pmFPAfileFringeIO.c \
    2526        pmFPAFlags.c \
    2627        pmFPALevel.c \
    … …  
    5152        pmFPAfileIO.h \
    5253        pmFPAfileFitsIO.h \
     54        pmFPAfileFringeIO.h \
    5355        pmFPAFlags.h \
    5456        pmFPALevel.h \

  • branches/pap/psModules/src/camera/pmFPABin.c

    r21183 r25027  
    2828    int numColsOut = binning->nXruff, numRowsOut = binning->nYruff; // Size of output image
    2929
    30     // Output image
    31     psImage *outImage;
     30
     31    psImage *outImage;                  // Output image
    3232    if (out->image && out->image->numCols >= numColsOut && out->image->numRows >= numRowsOut) {
    3333        outImage = out->image;
    3434    } else {
    3535        outImage = out->image = psImageRecycle(out->image,  numColsOut, numRowsOut, PS_TYPE_F32);
     36    }
     37
     38    psImage *outMask;                   // Output mask
     39    if (out->mask && out->mask->numCols >= numColsOut && out->mask->numRows >= numRowsOut) {
     40        outMask = out->mask;
     41    } else {
     42        outMask = out->mask = psImageRecycle(out->mask,  numColsOut, numRowsOut, PS_TYPE_IMAGE_MASK);
    3643    }
    3744
    … …  
    5360                        continue;
    5461                    }
     62                    if (!isfinite(inImage->data.F32[y][x])) {
     63                        continue;
     64                    }
    5565                    sum += inImage->data.F32[y][x];
    5666                    numPix++;
    … …  
    5868            }
    5969
    60             outImage->data.F32[yOut][xOut] = numPix > 0 ? sum / numPix : NAN;
     70            // Values to set
     71            float imageValue;
     72            psImageMaskType maskValue;
     73            if (numPix > 0) {
     74                imageValue = sum / numPix;
     75                maskValue = 0;
     76            } else {
     77                imageValue = NAN;
     78                maskValue = maskVal;
     79            }
     80            outImage->data.F32[yOut][xOut] = imageValue;
     81            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
    6182            xStart = xStop;
    6283        }

  • branches/pap/psModules/src/camera/pmFPACopy.c

    r23761 r25027  
    7070    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
    7171    *target = psImageAlloc(binning->nXruff, binning->nYruff, source->type.type);
     72    psImageInit (*target, 0.0);
    7273    return;
    7374}
    … …  
    221222    psTrace("psModules.camera", 3, "xFlip: %d; yFlip: %d\n", xFlip, yFlip);
    222223
     224    // Blow away extant readouts
     225    for (int i = 0; i < target->readouts->n; i++) {
     226        psFree(target->readouts->data[i]);
     227        target->readouts->data[i] = NULL;
     228    }
     229    target->readouts->n = 0;
     230
    223231    // Perform deep copy of the images.  I would prefer *not* to do a deep copy, in the interests of speed (we
    224232    // still need to do another deep copy into the HDU for when we write out), but this is the only way I can
    … …  
    242250        readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip,
    243251                             pixels);
     252        // Copy covariance matrix: doesn't care about flips, etc.
     253        if (sourceReadout->covariance) {
     254            if (targetReadout->covariance) {
     255                psFree(targetReadout->covariance);
     256            }
     257            targetReadout->covariance = psKernelCopy(sourceReadout->covariance);
     258#if 0
     259            if (binning) {
     260                // XXX This isn't strictly correct, but we don't have a function that bins covariance matrices
     261                // with unequal binning factors.
     262                psKernel *covar = psImageCovarianceBin(PS_MAX(binning->nXbin, binning->nYbin),
     263                                                       targetReadout->covariance);
     264                psFree(targetReadout->covariance);
     265                targetReadout->covariance = covar;
     266            }
     267#endif
     268        }
    244269
    245270        // Copy bias

  • branches/pap/psModules/src/camera/pmFPAMaskWeight.c

    r23259 r25027  
    109109    float saturation = psMetadataLookupF32(&mdok, cell->concepts, "CELL.SATURATION"); // Saturation level
    110110    if (!mdok || isnan(saturation)) {
    111         psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n");
    112         return false;
     111        // psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n");
     112        // return false;
     113        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
     114        saturation = NAN;
    113115    }
    114116    float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
    115117    if (!mdok || isnan(bad)) {
    116         psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n");
    117         return false;
     118        // psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n");
     119        // return false;
     120        psWarning("CELL.BAD is not set --- completely masking cell.\n");
     121        bad = NAN;
    118122    }
    119123    psTrace("psModules.camera", 5, "Saturation: %f, bad: %f\n", saturation, bad);
    120124
     125    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN;
     126    // in this case, we have to set the mask as well
     127    float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
     128    if (!mdok) { gain = NAN; }
     129    float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
     130    if (!mdok) { readnoise = NAN; }
    121131
    122132    // Set up the mask
    … …  
    127137    }
    128138    psImage *mask = readout->mask;      // The mask pixels
     139
     140    // completely mask if SATURATION or BAD are invalid
     141    if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) {
     142        psImageInit(mask, badMask);
     143        return true;
     144    }
     145
    129146    psImageInit(mask, 0);
    130147
    … …  
    199216}
    200217
    201 bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
     218bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson)
    202219{
    203220    PS_ASSERT_PTR_NON_NULL(readout, false);
     221    // check that the noiseMap (if it exists) matches the readout variance size)
    204222
    205223    pmCell *cell = readout->parent;     // The parent cell
    … …  
    209227    float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
    210228    if (!mdok || isnan(gain)) {
    211         psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n");
    212         return false;
     229        // psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n");
     230        // return false;
     231        psWarning("CELL.GAIN is not set --- setting variance to NAN\n");
     232        gain = NAN;
    213233    }
    214234    float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
    215235    if (!mdok || isnan(readnoise)) {
    216         psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n");
    217         return false;
    218     }
    219     if (psMetadataLookup(cell->concepts, "CELL.READNOISE.UPDATE")) {
     236        // psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n");
     237        // return false;
     238        psWarning("CELL.READNOISE is not set --- setting variance to NAN\n");
     239        readnoise = NAN;
     240    }
     241    // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN)
     242    if (!isnan(gain) && psMetadataLookup(cell->concepts, "CELL.READNOISE.UPDATE")) {
    220243        psError(PS_ERR_IO, true, "CELL.READNOISE has not yet been updated for the gain");
    221244        return false;
     245    }
     246
     247    // for invalid input data, set the readout variance to NAN
     248    if (isnan(gain) || isnan(readnoise)) {
     249        if (!readout->variance) {
     250            // generate the image if needed
     251            readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
     252        }
     253        // XXX need to set the mask, if defined
     254        psImageInit(readout->variance, NAN);
     255        return true;
    222256    }
    223257
    … …  
    228262
    229263        // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
     264        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
    230265        readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
    231266        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
    … …  
    239274    }
    240275
    241     readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+",
    242                                            psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     276    // apply a supplied readnoise map (NOTE: in DN, not electrons):
     277    if (noiseMap) {
     278        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
     279        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
     280        psFree (rdVar);
     281    } else {
     282        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     283    }
    243284
    244285    return true;
    … …  
    247288// this function creates the variance pixels, or uses the existing variance pixels.  it will set
    248289// the noise pixel values only if the variance image is not supplied
    249 bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson)
     290bool pmReadoutGenerateVariance(pmReadout *readout, const psImage *noiseMap, bool poisson)
    250291{
    251292    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    291332    readout->variance = variance;
    292333
    293     return pmReadoutSetVariance(readout, poisson);
    294 }
    295 
    296 bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     334    return pmReadoutSetVariance(readout, noiseMap, poisson);
     335}
     336
     337bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson)
    297338{
    298339    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    301342
    302343    success &= pmReadoutGenerateMask(readout, satMask, badMask);
    303     success &= pmReadoutGenerateVariance(readout, poisson);
     344    success &= pmReadoutGenerateVariance(readout, noiseMap, poisson);
    304345
    305346    return success;
    306347}
    307348
    308 bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     349bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson)
    309350{
    310351    PS_ASSERT_PTR_NON_NULL(cell, false);
    … …  
    314355    for (int i = 0; i < readouts->n; i++) {
    315356        pmReadout *readout = readouts->data[i]; // The readout
    316         success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask);
     357        success &= pmReadoutGenerateMaskVariance(readout, satMask, badMask, noiseMap, poisson);
    317358    }
    318359
    … …  
    482523        return false;
    483524    }
    484     float stdev = psStatsGetValue(stdevStats, stdevStat); // Stadard deviation of fluxes
     525    float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes
    485526    psFree(stdevStats);
    486527    psFree(noise);

  • branches/pap/psModules/src/camera/pmFPAMaskWeight.h

    r21363 r25027  
    5454/// can't be generated.
    5555bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance
     56                          const psImage *noiseMap, ///< 2D image of the read noise in DN
    5657                          bool poisson    ///< Include poisson variance (in addition to read noise)?
    5758    );
    … …  
    7273/// with HDU entry).  This is intended for most operations.
    7374bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance
     75                          const psImage *noiseMap, ///< 2D image of the read noise in DN
    7476                               bool poisson    ///< Include poisson variance (in addition to read noise)?
    7577    );
    … …  
    8183                                   psImageMaskType sat, ///< Mask value to give saturated pixels
    8284                                   psImageMaskType bad, ///< Mask value to give bad (low) pixels
     85                                   const psImage *noiseMap, ///< 2D image of the read noise in DN
    8386                                   bool poisson ///< Include poisson variance (in addition to read noise)?
    8487    );
    … …  
    9093                                psImageMaskType sat, ///< Mask value to give saturated pixels
    9194                                psImageMaskType bad, ///< Mask value to give bad (low) pixels
     95                                const psImage *noiseMap, ///< 2D image of the read noise in DN
    9296                                bool poisson ///< Include poisson variance (in addition to read noise)?
    9397    );

  • branches/pap/psModules/src/camera/pmFPARead.c

    r23827 r25027  
    6464
    6565// Set the "thisXXXScan" value in the readout for the appropriate image type
    66 static int readoutSetThisScan(pmReadout *readout, // Readout of interest
     66static void readoutSetThisScan(pmReadout *readout, // Readout of interest
    6767                              fpaReadType type, // Type of image
    6868                              int thisScan // Starting scan number
    … …  
    7272      case FPA_READ_TYPE_IMAGE:
    7373        readout->thisImageScan = thisScan;
    74         return readout->lastImageScan;
     74        return;
    7575      case FPA_READ_TYPE_MASK:
    7676        readout->thisMaskScan = thisScan;
    77         return readout->lastMaskScan;
     77        return;
    7878      case FPA_READ_TYPE_VARIANCE:
    7979        readout->thisVarianceScan = thisScan;
    80         return readout->lastVarianceScan;
     80        return;
    8181      default:
    8282        psAbort("Unknown read type: %x\n", type);
    8383    }
    84     return false;
    8584}
    8685
    … …  
    103102
    104103// Set the "lastXXXScan" value in the readout for the appropriate image type
    105 static int readoutSetLastScan(pmReadout *readout, // Readout of interest
     104static void readoutSetLastScan(pmReadout *readout, // Readout of interest
    106105                              fpaReadType type, // Type of image
    107106                              int lastScan // Last scan number
    … …  
    111110      case FPA_READ_TYPE_IMAGE:
    112111        readout->lastImageScan = lastScan;
    113         return readout->lastImageScan;
     112        return;
    114113      case FPA_READ_TYPE_MASK:
    115114        readout->lastMaskScan = lastScan;
    116         return readout->lastMaskScan;
     115        return;
    117116      case FPA_READ_TYPE_VARIANCE:
    118117        readout->lastVarianceScan = lastScan;
    119         return readout->lastVarianceScan;
     118        return;
    120119      default:
    121120        psAbort("Unknown read type: %x\n", type);
    122121    }
    123     return false;
    124122}
    125123
    … …  
    143141// Determine number of readouts in the FITS file
    144142// In the process, reads the header and concepts
    145 static bool cellNumReadouts(pmCell *cell,    // Cell of interest
     143static int cellNumReadouts(pmCell *cell,    // Cell of interest
    146144                            psFits *fits,    // FITS file
    147145                            pmConfig *config // Configuration
    … …  
    155153    if (!hdu || hdu->blankPHU) {
    156154        psError(PS_ERR_IO, true, "Unable to find HDU");
    157         return false;
     155        return 0;
    158156    }
    159157    if (!pmCellReadHeader(cell, fits, config)) {
    160158        psError(PS_ERR_IO, false, "Unable to read header for cell!\n");
    161         return false;
     159        return 0;
    162160    }
    163161    if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS, true, NULL)) {
    164162        psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n");
    165         return false;
     163        return 0;
    166164    }
    167165
    … …  
    171169    if (!mdok) {
    172170        psError(PS_ERR_IO, true, "Unable to find NAXIS in header for extension %s\n", hdu->extname);
    173         return false;
    174     }
     171        return 0;
     172    }
     173
    175174    if (naxis == 0) {
    176175        // No pixels to read
    177176        psError(PS_ERR_IO, true, "No pixels in extension %s.", hdu->extname);
    178         return false;
     177        return 0;
    179178    }
    180179    if (naxis < 2 || naxis > 3) {
    181180        psError(PS_ERR_IO, true, "NAXIS in header of extension %s (= %d) is not valid.\n",
    182181                hdu->extname, naxis);
    183         return false;
     182        return 0;
    184183    }
    185184    int naxis3;                     // Number of image planes
    … …  
    188187        if (!mdok) {
    189188            psError(PS_ERR_IO, true, "Unable to find NAXIS3 in header for extension %s\n", hdu->extname);
    190             return false;
     189            return 0;
    191190        }
    192191    } else {
    … …  
    296295static bool readoutMore(pmReadout *readout, // Readout of interest
    297296                        psFits *fits,    // FITS file
    298                         int z,          // Plane number
     297                        int z,          // Plane number to read
     298                        int *zMax,      // Max plane number in this cell
    299299                        int numScans,   // Number of scans to read at a time
    300300                        fpaReadType type, // Type of image
    … …  
    310310    // N readouts, but numScans set to 0.  only the first should report that it requires data,
    311311    // even if all readouts lack the image pointer.
    312     if (!image) {
     312    if (numScans == 0) {
     313      if (!image) {
    313314        return true;
    314     }
    315 
    316     // If we have already read an image, this result implies we are done (no more to read)
    317     if (numScans == 0) {
    318         return false;
     315      } else {
     316        return false;
     317      }
    319318    }
    320319
    … …  
    324323        return false;
    325324    }
    326     int naxis3 = cellNumReadouts(cell, fits, config); // Number of planes
    327     if (z >= naxis3) {
     325    *zMax = cellNumReadouts(cell, fits, config); // Number of planes
     326    if (z >= *zMax) {
    328327        // No more to read
    329328        return false;
    … …  
    486485                             psFits *fits, // FITS file
    487486                             int z,     // Desired image plane
     487                             int *zMax, // Max plane number in this cell
    488488                             int numScans, // Number of scans (row or col depends on CELL.READDIR); 0 for all
    489489                             int overlap, // Number of scans (row/col) to overlap between scans
    … …  
    511511
    512512    int naxis3 = cellNumReadouts(cell, fits, config); // Number of image planes
     513    if (zMax) *zMax = naxis3;
     514
    513515    if (z >= naxis3) {
    514516        psError(PS_ERR_IO, false, "Desired image plane (%d) exceeds available number (%d).",
    … …  
    592594    }
    593595
    594     // Determine the number of scans to read
    595     int lastScan = readoutSetLastScan(readout, type, thisScan + numScans);
     596    int origThisScan = thisScan;        // Original value of thisScan (starting point for read)
    596597    if (thisScan == 0) {
    597598        overlap = 0;
    … …  
    601602        thisScan = 0;
    602603    }
    603     readoutSetThisScan(readout, type, thisScan);
    604604
    605605    // Calculate limits, adjust readout->row0,col0
    … …  
    620620        }
    621621    }
     622    int lastScan = origThisScan + numScans; // Last scan to read
     623
     624    readoutSetThisScan(readout, type, thisScan);
     625    readoutSetLastScan(readout, type, lastScan);
    622626
    623627    // Blow away existing data.
    … …  
    10211025
    10221026
    1023 bool pmReadoutMore(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
     1027bool pmReadoutMore(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config)
    10241028{
    10251029    PS_ASSERT_PTR_NON_NULL(readout, false);
    10261030    PS_ASSERT_FITS_NON_NULL(fits, false);
    10271031
    1028     return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_IMAGE, config);
    1029 }
    1030 
    1031 bool pmReadoutReadChunk(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, pmConfig *config)
     1032    return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_IMAGE, config);
     1033}
     1034
     1035bool pmReadoutReadChunk(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap, pmConfig *config)
    10321036{
    10331037    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    10361040    PS_ASSERT_INT_NONNEGATIVE(numScans, false);
    10371041
    1038     return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_IMAGE, config);
     1042    return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_IMAGE, config);
    10391043}
    10401044
    … …  
    10441048    PS_ASSERT_FITS_NON_NULL(fits, false);
    10451049
    1046     return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_IMAGE, config);
     1050    return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_IMAGE, config);
    10471051}
    10481052
    … …  
    10841088//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    10851089
    1086 bool pmReadoutMoreMask(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
     1090bool pmReadoutMoreMask(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config)
    10871091{
    10881092    PS_ASSERT_PTR_NON_NULL(readout, false);
    10891093    PS_ASSERT_FITS_NON_NULL(fits, false);
    10901094
    1091     return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_MASK, config);
    1092 }
    1093 
    1094 bool pmReadoutReadChunkMask(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
     1095    return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_MASK, config);
     1096}
     1097
     1098bool pmReadoutReadChunkMask(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap,
    10951099                            pmConfig *config)
    10961100{
    … …  
    11001104    PS_ASSERT_INT_NONNEGATIVE(numScans, false);
    11011105
    1102     return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_MASK, config);
     1106    return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_MASK, config);
    11031107}
    11041108
    … …  
    11081112    PS_ASSERT_FITS_NON_NULL(fits, false);
    11091113
    1110     return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_MASK, config);
     1114    return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_MASK, config);
    11111115}
    11121116
    … …  
    11391143//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    11401144
    1141 bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
     1145bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config)
    11421146{
    11431147    PS_ASSERT_PTR_NON_NULL(readout, false);
    11441148    PS_ASSERT_FITS_NON_NULL(fits, false);
    11451149
    1146     return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config);
    1147 }
    1148 
    1149 bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
     1150    return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_VARIANCE, config);
     1151}
     1152
     1153bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap,
    11501154                              pmConfig *config)
    11511155{
    … …  
    11551159    PS_ASSERT_INT_NONNEGATIVE(numScans, false);
    11561160
    1157     return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config);
     1161    return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_VARIANCE, config);
    11581162}
    11591163
    … …  
    11631167    PS_ASSERT_FITS_NON_NULL(fits, false);
    11641168
    1165     return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config);
     1169    return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_VARIANCE, config);
    11661170}
    11671171

  • branches/pap/psModules/src/camera/pmFPARead.h

    r21363 r25027  
    2121                   psFits *fits,        ///< FITS file from which to read
    2222                   int z,               ///< Readout number/plane; zero-offset indexing
     23                   int *zMax,           ///< Max plane number in this cell
    2324                   int numScans,        ///< Number of scans (rows/cols) to read
    2425                   pmConfig *config     ///< Configuration
    … …  
    3132                        psFits *fits,   ///< FITS file from which to read
    3233                        int z,          ///< Readout number/plane; zero-offset indexing
     34                   int *zMax,           ///< Max plane number in this cell
    3335                        int numScans,   ///< Number of scans (rows/cols) to read
    3436                        int overlap,    ///< Overlap between consecutive reads
    … …  
    98100                       psFits *fits,    ///< FITS file from which to read
    99101                       int z,           ///< Readout number/plane; zero-offset indexing
     102                       int *zMax,       ///< Max plane number in this cell
    100103                       int numScans,    ///< Number of scans (rows/cols) to read
    101104                       pmConfig *config ///< Configuration
    … …  
    108111                            psFits *fits, ///< FITS file from which to read
    109112                            int z,      ///< Readout number/plane; zero-offset indexing
     113                            int *zMax,          ///< Max plane number in this cell
    110114                            int numScans, ///< Number of scans (rows/cols) to read
    111115                            int overlap, ///< Overlap between consecutive reads
    … …  
    150154                           psFits *fits, ///< FITS file from which to read
    151155                           int z,       ///< Readout number/plane; zero-offset indexing
     156                           int *zMax,   ///< Max plane number in this cell
    152157                           int numScans, ///< Number of scans (rows/cols) to read
    153158                           pmConfig *config ///< Configuration
    … …  
    160165                                psFits *fits, ///< FITS file from which to read
    161166                                int z,    ///< Readout number/plane; zero-offset indexing
     167                   int *zMax,           ///< Max plane number in this cell
    162168                                int numScans, ///< Number of scans (rows/cols) to read
    163169                                int overlap, ///< Overlap between consecutive reads

  • branches/pap/psModules/src/camera/pmFPAWrite.c

    r23428 r25027  
    160160    psArray **imageArray = appropriateImageArray(hdu, type); // Array of images in the HDU
    161161
     162    // XXX detect missing variance & mask images...
     163
    162164    // Generate the HDU if needed --- this is required after a pmFPACopy, or similar, which does not
    163165    // generate the HDU, but only copies the structure.
    164     if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) {
    165         psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.");
    166         return false;
     166    if (!blank && !hdu->blankPHU && !*imageArray) {
     167        if (!pmHDUGenerateForCell(cell)) {
     168            psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.");
     169            return false;
     170        }
     171        if (!*imageArray) {
     172            if (type == FPA_WRITE_TYPE_IMAGE) {
     173                psError(PS_ERR_UNKNOWN, false, "Expected to write an image, but it is missing...programming error?.");
     174                return false;
     175            }
     176            if (type == FPA_WRITE_TYPE_MASK) {
     177                psWarning("No mask image for this cell; skipping");
     178            }
     179            if (type == FPA_WRITE_TYPE_VARIANCE) {
     180                psWarning("No variance image for this cell; skipping");
     181            }
     182            return true;
     183        }
    167184    }
    168185
    … …  
    724741    }
    725742    if (!psFitsCompressionApply(fits, compress)) {
    726         psError(PS_ERR_UNKNOWN, false, "Unable to set FITS compression to NONE");
     743        psError(PS_ERR_UNKNOWN, false, "Unable to restore FITS compression");
    727744        psFree(compress);
    728745        return false;

  • branches/pap/psModules/src/camera/pmFPAfileIO.c

    r23576 r25027  
    2323#include "pmFPAWrite.h"
    2424#include "pmFPAfileFitsIO.h"
     25#include "pmFPAfileFringeIO.h"
    2526#include "pmSpan.h"
    2627#include "pmFootprint.h"
    … …  
    192193        status = pmFPAviewReadFitsImage(view, file, config);
    193194        if (status) {
    194             if (!pmFPAviewReadFitsTable(view, file, "FRINGE")) {
     195            if (!pmFPAviewReadFringes(view, file)) {
    195196                psError(PS_ERR_UNKNOWN, false, "Unable to read fringe data from %s.\n", file->filename);
    196197                return false;
    … …  
    451452        status = pmFPAviewWriteFitsImage (view, file, config);
    452453        if (status) {
    453             if (!pmFPAviewWriteFitsTable(view, file, "FRINGE", config)) {
     454            if (!pmFPAviewWriteFringes(view, file, config)) {
    454455                psError(PS_ERR_UNKNOWN, false, "Unable to write fringe data from %s.\n", file->filename);
    455456                return false;

  • branches/pap/psModules/src/camera/pmHDU.c

    r21363 r25027  
    9393    }
    9494
     95    psTrace("psModules.camera", 5, "Reading the header...\n");
     96
     97    // The header may already exist (e.g., from doing concept writing at the PHU level) so we need to be
     98    // careful.  We read into a separate container and copy that over the top of anything that's already read.
     99    psMetadata *header = psFitsReadHeader(NULL, fits);
     100    if (!header) {
     101        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname);
     102        return false;
     103    }
     104
    95105    if (!hdu->header) {
    96         psTrace("psModules.camera", 5, "Reading the header...\n");
    97         hdu->header = psFitsReadHeader(hdu->header, fits);
    98         if (! hdu->header) {
    99             psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname);
    100             return false;
    101         }
    102     }
     106        hdu->header = header;
     107        return true;
     108    }
     109
     110    psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator
     111    psMetadataItem *item;           // Item from iteration
     112    while ((item = psMetadataGetAndIncrement(iter))) {
     113        const char *name = item->name; // Name of item
     114        if (psMetadataLookup(hdu->header, name)) {
     115            // It exists; clobber
     116            psMetadataRemoveKey(hdu->header, name);
     117        }
     118        psMetadataAddItem(hdu->header, item, PS_LIST_TAIL, 0);
     119    }
     120    psFree(iter);
     121    psFree(header);
    103122
    104123    return true;

  • branches/pap/psModules/src/camera/pmHDUGenerate.c

    r23318 r25027  
    390390    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) {
    391391        // Nothing from which to create an HDU
    392         psFree(cells);
    393         psError(PS_ERR_IO, true, "Nothing from which to create an HDU\n");
    394         return false;
     392        // psError(PS_ERR_IO, true, "Nothing from which to create an HDU\n");
     393        psWarning("Nothing from which to create an HDU, must be empty\n");
     394        return true;
    395395    }
    396396
    … …  
    504504        psFree(iter);
    505505    }
    506     psFree(cells);
    507 
    508506    return true;
    509507}
    … …  
    615613                return true;
    616614            }
    617 
    618             return generateHDU(hdu, cells);
     615            bool status = generateHDU(hdu, cells);
     616            psFree (cells);
     617            return status;
    619618        }
    620619    case PM_FPA_LEVEL_CHIP:
    … …  
    664663            }
    665664
    666             return generateHDU(hdu, cells);
     665            bool status = generateHDU(hdu, cells);
     666            psFree (cells);
     667            return status;
    667668        }
    668669    case PM_FPA_LEVEL_FPA:

  • branches/pap/psModules/src/camera/pmReadoutFake.c

    r21104 r25027  
    2929#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
    3030#define SOURCE_MASK (PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
     31#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
     32                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
    3133
    3234
    … …  
    8082
    8183    int numSources = sources->n;          // Number of stars
    82 
    83     float flux0 = NAN;                  // Flux for central intensity of 1.0
    84     if (normalisePeak) {
    85         pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0,
    86                                                  1.0); // Fake model, with central intensity of 1.0
    87         psAssert(fakeModel, "failed to generate model: should this be an error or not?");
    88 
    89         if (circularise && !circulariseModel(fakeModel)) {
    90             psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
    91             psFree(fakeModel);
    92             return false;
    93         }
    94 
    95         flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0
    96         psFree(fakeModel);
    97     }
    98 
    9984    for (int i = 0; i < numSources; i++) {
    10085        pmSource *source = sources->data[i]; // Source of interest
    … …  
    118103
    119104        float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source
     105
    120106        if (normalisePeak) {
    121             flux /= flux0;
     107            // Normalise flux
     108            pmModel *normModel = pmModelFromPSFforXY(psf, x, y, 1.0); // Model for normalisation
     109            if (!normModel || (normModel->flags & MODEL_MASK)) {
     110                psFree(normModel);
     111                continue;
     112            }
     113            if (circularise && !circulariseModel(normModel)) {
     114                psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
     115                psFree(normModel);
     116                return false;
     117            }
     118
     119            flux /= normModel->modelFlux(normModel->params);
     120            psFree(normModel);
    122121        }
    123122
    124123        pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, flux);
    125         if (!fakeModel) {
     124        if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
     125            psFree(fakeModel);
    126126            continue;
    127127        }

  • branches/pap/psModules/src/camera/pmReadoutStack.c

    r21363 r25027  
    252252            continue;
    253253        }
     254        if (!readout->process) {
     255            continue;
     256        }
    254257        if (!readout->image) {
    255258            psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %ld has NULL image.\n", i);
    … …  
    260263        pmCell *cell = readout->parent; // The parent cell
    261264        bool mdok = true;       // Status of MD lookup
    262         psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
     265        psRegion *trimsec = cell ? psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC") : NULL; // Trim section
    263266        if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) {
    264             psWarning("CELL.TRIMSEC is not set for readout %ld --- ignored.\n", i);
     267            psWarning("CELL.TRIMSEC is not set for readout %ld --- attempting to use image size.\n", i);
     268            xSize = PS_MAX(xSize, readout->image->numCols);
     269            ySize = PS_MAX(ySize, readout->image->numRows);
     270            xMin = PS_MIN(xMin, 0);
     271            xMax = PS_MAX(xMax, readout->image->numCols - 1);
     272            yMin = PS_MIN(yMin, 0);
     273            yMax = PS_MAX(yMax, readout->image->numRows - 1);
    265274        } else {
    266275            xSize = PS_MAX(xSize, trimsec->x1 - trimsec->x0);

  • branches/pap/psModules/src/concepts/pmConcepts.c

    r23832 r25027  
    485485        concept[length - 1] = '\0';
    486486
     487        // special variants:
     488        if (!strcmp(concept, "FPA.DATE")) {
     489          psTime *fpaTime = psMetadataLookupPtr(NULL, fpa->concepts, "FPA.TIME");
     490          if (!fpaTime) {
     491            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Missing concept FPA.TIME needed for FPA.DATE");
     492            psFree(string);
     493            return NULL;
     494          }
     495          psString dateTimeString = psTimeToISO(fpaTime); // String representation
     496          psList *dateTime = psStringSplit(dateTimeString, "T", true);
     497          psFree(dateTimeString);
     498          psString dateString = psMemIncrRefCounter(psListGet(dateTime, PS_LIST_HEAD)); // The date string
     499          psFree (dateTime);
     500
     501          if (!psStringSubstitute(&string, dateString, "{FPA.DATE}")) {
     502              psError(PS_ERR_UNKNOWN, false, "Unable to replace concept %s", concept);
     503              psFree(string);
     504              psFree(dateString);
     505              return NULL;
     506          }
     507          psFree (dateString);
     508          continue;
     509        }
     510
    487511        psTrace("psModules.concepts", 7, "Interpolating concept %s", concept);
    488512

  • branches/pap/psModules/src/concepts/pmConceptsRead.c

    r22699 r25027  
    425425    if (rows->n == 0) {
    426426        psWarning("No rows returned from database query for concept %s --- ignored.", name);
     427        psFree(rows);
    427428        return NULL;
    428429    }

  • branches/pap/psModules/src/concepts/pmConceptsStandard.c

    r23819 r25027  
    3636// Format type for time
    3737typedef enum {
    38   TIME_FORMAT_YYYYMMDD,                 // Date stored in YYYY-MM-DD order (ISO-standard)
    39   TIME_FORMAT_DDMMYYYY,                 // Date stored in DD-MM-YYYY order
    40   TIME_FORMAT_MMDDYYYY,                 // Date stored in MM-DD-YYYY order
    41   TIME_FORMAT_JD,                       // Date stored as JD
    42   TIME_FORMAT_MJD,                      // Date stored as MJD
     38  TIME_FORMAT_YYYYMMDD,                 // Date stored in YYYY-MM-DD order (ISO-standard)
     39  TIME_FORMAT_DDMMYYYY,                 // Date stored in DD-MM-YYYY order
     40  TIME_FORMAT_MMDDYYYY,                 // Date stored in MM-DD-YYYY order
     41  TIME_FORMAT_JD,                       // Date stored as JD
     42  TIME_FORMAT_MJD,                      // Date stored as MJD
    4343} conceptTimeFormatType;
    4444
    … …  
    391391}
    392392
     393
    393394// FPA.RA and FPA.DEC
    394395psMetadataItem *p_pmConceptFormat_FPA_Coords(const psMetadataItem *concept,
    … …  
    410411    // How to interpret the coordinates
    411412    bool mdok = true;                   // Status of MD lookup
    412     psMetadata *formats = psMetadataLookupMetadata(&mdok,
    413                                                    cameraFormat,
    414                                                    "FORMATS");
     413    psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS");
     414    bool sexagesimal = false;           // Write sexagesimal format?
    415415    if (mdok && formats) {
    416         psString format = psMetadataLookupStr(&mdok,formats, concept->name);
     416        psString format = psMetadataLookupStr(&mdok, formats, concept->name);
    417417        if (mdok && strlen(format) > 0) {
    418418            if (strcasecmp(format, "HOURS") == 0) {
    … …  
    428428            coords /= defaultCoordScaling(concept);
    429429        }
     430
     431        psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA
     432        psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec
     433        if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) {
     434            sexagesimal = true;
     435        }
    430436    } else {
    431437        coords /= defaultCoordScaling(concept);
    432438    }
    433439
    434     // We choose to write sexagesimal format
    435     int big, medium;                    // Degrees and minutes
    436     float small;                        // Seconds
    437     bool negative = (coords < 0);       // Are we working below zero?
    438     coords = fabs(coords);
    439     big = (int)abs(coords);
    440     coords -= big;
    441     medium = 60.0 * coords;
    442     coords -= medium / 60.0;
    443     small = 3600.0 * coords;
    444     small = (float)((int)(small * 1000.0)) / 1000.0;
    445     if (negative) {
    446         big *= -1;
    447     }
    448     psString coordString = NULL;        // String with the coordinates in sexagesimal format
    449     psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small);
    450     psMetadataItem *coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString);
    451     psFree(coordString);
     440    psMetadataItem *coordItem = NULL;   // Item with coordinates, to return
     441    if (sexagesimal) {
     442        int big, medium;                    // Degrees and minutes
     443        float small;                        // Seconds
     444        bool negative = (coords < 0);       // Are we working below zero?
     445        coords = fabs(coords);
     446        big = (int)abs(coords);
     447        coords -= big;
     448        medium = 60.0 * coords;
     449        coords -= medium / 60.0;
     450        small = 3600.0 * coords;
     451        small = (float)((int)(small * 1000.0)) / 1000.0;
     452        if (negative) {
     453            big *= -1;
     454        }
     455        psString coordString = NULL;        // String with the coordinates in sexagesimal format
     456        psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small);
     457        coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString);
     458        psFree(coordString);
     459    } else {
     460        coordItem = psMetadataItemAllocF64(concept->name, concept->comment, coords);
     461    }
    452462
    453463    return coordItem;

  • branches/pap/psModules/src/config/pmConfig.c

    r23748 r25027  
    243243    unsigned int numBadLines = 0;
    244244    struct stat filestat;
     245
     246    pmErrorRegister();
    245247
    246248    psTrace("psModules.config", 3, "Loading %s configuration from file %s\n",
    … …  
    14961498                }
    14971499            }
     1500            globfree(&globList);
    14981501        }
    14991502        psFree (words);

  • branches/pap/psModules/src/config/pmErrorCodes.dat

    r23688 r25027  
    1212SKY                     Problem in sky
    1313STAMPS                  Unable to select stamps for PSF-matching
     14SMALL_AREA              Small area for PSF-matching
    1415# these errors correspond to standard exit conditions
    1516ARGUMENTS               Incorrect arguments

  • branches/pap/psModules/src/detrend/Makefile.am

    r20617 r25027  
    1616        pmShifts.c \
    1717        pmDark.c \
    18         pmRemnance.c
     18        pmRemnance.c \
     19        pmPattern.c
    1920
    2021#       pmSkySubtract.c
    … …  
    3334        pmShifts.h \
    3435        pmDark.h \
    35         pmRemnance.h
     36        pmRemnance.h \
     37        pmPattern.h
    3638
    3739#       pmSkySubtract.h

  • branches/pap/psModules/src/detrend/pmDark.c

    r21183 r25027  
     1#ifdef HAVE_CONFIG_H
     2#include <config.h>
     3#endif
     4
    15#include <stdio.h>
    26#include <pslib.h>
    37#include <string.h>
     8#include <strings.h>
    49
    510#include "psPolynomialMD.h"
    … …  
    1217#include "pmReadoutStack.h"
    1318#include "pmDetrendThreads.h"
    14 
     19#include "pmErrorCodes.h"
    1520#include "pmDark.h"
    1621
    1722#define PM_DARK_FITS_EXTNAME "PS_DARK"  // FITS extension name for ordinates table
    1823#define PM_DARK_FITS_NAME    "NAME"     // Column name for concept name in ordinates table
     24#define PM_DARK_FITS_RULE    "RULE"     // Column name for concept rule in ordinates table
    1925#define PM_DARK_FITS_ORDER   "ORDER"    // Column name for polynomial order in ordinates table
    2026#define PM_DARK_FITS_SCALE   "SCALE"    // Column name for scaling option in ordinates table
    … …  
    2228#define PM_DARK_FITS_MAX     "MAX"      // Column name for maximum value in ordinates table
    2329
    24 
    25 
    26 // Look up the value of an ordinate in a readout
    27 static bool ordinateLookup(float *value, // Value of ordinate, to return
    28                            bool *inRange, // is value within min : max range?
    29                            const char *name, // Name of ordinate (concept name)
    30                            bool scale,  // Scale the value?
    31                            float min, float max, // Minimum and maximum values for scaling
    32                            const pmReadout *ro // Readout of interest
    33                            )
    34 {
    35     assert(value);
    36     assert(name);
    37     assert(ro);
    38 
    39     *inRange = true;
    40 
    41     pmCell *cell = ro->parent; // Parent cell
    42     if (!cell) {
    43         return false;
    44     }
     30static bool ordinateParseConcept(double *value, const pmReadout *readout, const char *name) {
     31
     32    *value = NAN;
     33
     34    pmCell *cell = readout->parent; // Parent cell
     35    psAssert(cell, "readout is missing cell \n");
     36
    4537    psMetadataItem *item = psMetadataLookup(cell->concepts, name);
    4638    if (!item) {
    4739        pmChip *chip = cell->parent; // Parent chip
    48         if (!chip) {
    49             return false;
    50         }
     40        psAssert(chip, "cell is missing chip \n");
     41
    5142        item = psMetadataLookup(chip->concepts, name);
    5243        if (!item) {
    5344            pmFPA *fpa = chip->parent; // Parent FPA
    54             if (!fpa) {
    55                 return false;
    56             }
     45            psAssert(fpa, "chip is missing fpa \n");
     46
    5747            item = psMetadataLookup(fpa->concepts, name);
    5848            if (!item) {
    59                 psWarning("Unable to find concept %s in readout", name);
     49                psError(PM_ERR_CONFIG, true, "Unable to find concept %s in readout", name);
    6050                return false;
    6151            }
    … …  
    6353    }
    6454
    65     *value = psMetadataItemParseF32(item); // Value of interest
     55    *value = psMetadataItemParseF64(item); // Value of interest
    6656    if (!isfinite(*value)) {
    6757        psWarning("Non-finite value (%f) of concept %s in readout", *value, name);
    68         return false;
    69     }
     58    }
     59    return true;
     60}
     61
     62static bool ordinateParseRule(double *value, const pmReadout *readout, const char *name, const char *rule) {
     63
     64    psAssert(name, "ordinate name is not defined");
     65    psAssert(rule, "ordinate rule is not defined");
     66
     67    *value = NAN;
     68
     69    psArray *words = psStringSplitArray(rule, " ", false);
     70   
     71    // we should have a rule of the form (concept) OP (concept) OP (concept) ...
     72    // for now, the only allowed OP is * (eventually, we can steal code from opihi for a better
     73    // RPN parser).
     74
     75    if (words->n % 2 == 0) {
     76        psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule);
     77        psFree(words);
     78        return false;
     79    }
     80
     81    for (int i = 1; i < words->n; i+=2) {
     82        if (strcmp((char *)words->data[i], "*")) {
     83            psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule);
     84            psFree(words);
     85            return false;
     86        }
     87    }
     88   
     89    if (!ordinateParseConcept(value, readout, words->data[0])) {
     90        psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule);
     91        psFree(words);
     92        return false;
     93    }
     94
     95    double value2 = 0.0;
     96    for (int i = 2; i < words->n; i+=2) {
     97        if (!ordinateParseConcept(&value2, readout, words->data[i])) {
     98            psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule);
     99            psFree(words);
     100            return false;
     101        }
     102        *value *= value2;
     103    }
     104    psFree(words);
     105    return true;
     106}
     107
     108// Look up the value of an ordinate in a readout
     109static bool ordinateLookup(double *value, // Value of ordinate, to return
     110                           bool *inRange, // is value within min : max range?
     111                           const char *name, // Name of ordinate (concept or abstract name)
     112                           const char *rule, // Rule for generating the value (if NULL use name as concept)
     113                           bool scale,  // Scale the value?
     114                           float min, float max, // Minimum and maximum values for scaling
     115                           const pmReadout *readout // Readout of interest
     116                           )
     117{
     118    assert(value);
     119    assert(name);
     120    assert(readout);
     121
     122    *inRange = true;
     123
     124    if (rule) {
     125        if (!ordinateParseRule(value, readout, name, rule)) {
     126            psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name);
     127            return false;
     128        }
     129    } else {
     130        if (!ordinateParseConcept(value, readout, name)) {
     131            psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name);
     132            return false;
     133        }
     134    }
     135
    70136    if (scale) {
    71137        if (*value < min || *value > max) {
    … …  
    82148{
    83149    psFree(ord->name);
     150    psFree(ord->rule);
    84151    return;
    85152}
    … …  
    91158
    92159    ord->name = psStringCopy(name);
     160    ord->rule = NULL;
    93161    ord->order = order;
    94162    ord->scale = false;
    … …  
    118186
    119187        pmReadout *readout = inputs->data[i]; // Readout of interest
    120         float normValue;            // Normalisation value
    121         if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) {
    122             psWarning("Unable to find value of %s for readout %d", normConcept, i);
     188        double normValue;            // Normalisation value
     189        if (!ordinateLookup(&normValue, &inRange, normConcept, NULL, false, NAN, NAN, readout)) {
     190            psError(PM_ERR_CONFIG, false, "problem finding concept %s for DARK.NORM", normConcept);
     191            return false;
     192        }
     193        if (!isfinite(normValue)) {
     194            psWarning("Unable to find acceptable value of %s for readout %d", normConcept, i);
    123195            roMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    124196            norm->data.F32[i] = NAN;
    … …  
    140212        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information
    141213        if (ord->order <= 0) {
    142             psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order);
     214            psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order);
    143215            psFree(values);
    144216            psFree(roMask);
    … …  
    157229
    158230            pmReadout *readout = inputs->data[j]; // Readout of interest
    159             float value = NAN;          // Value of ordinate
    160             if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) {
     231            double value = NAN;          // Value of ordinate
     232            if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) {
     233                psError(PM_ERR_CONFIG, false, "problem finding rule for DARK.ORDINATE %s", ord->name);
     234                return false;
     235            }
     236            if (!isfinite(value)) {
     237                psWarning("Unable to find acceptable value of DARK.ORDINATE %s for readout %d", ord->name, i);
    161238                roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
    162239                val->data.F32[i] = NAN;
    … …  
    165242            }
    166243            if (!inRange) {
     244                psWarning("Value of DARK.ORDINATE %s for readout %d is out of range", ord->name, i);
    167245                roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
    168246                val->data.F32[i] = NAN;
    … …  
    309387        return false;
    310388    }
     389
     390    pmDarkVisualInit(values);
    311391
    312392    pmReadout *outReadout = output->readouts->data[0];
    … …  
    352432                psVectorInit(poly->coeff, NAN);
    353433            }
     434
     435            pmDarkVisualPixelFit(pixels, mask);
     436            pmDarkVisualPixelModel(poly, values);
     437
    354438            for (int k = 0; k < poly->coeff->n; k++) {
    355439                pmReadout *ro = output->readouts->data[k]; // Readout of interest
    … …  
    447531    for (int i = 0; i < numOrdinates; i++) {
    448532        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
    449         float value = NAN;              // Value for ordinate
    450         if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) {
    451             psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", ord->name);
     533        double value = NAN;              // Value for ordinate
     534        if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) {
     535            psError(PS_ERR_UNKNOWN, true, "Unable to find value for DARK.ORDINATE %s", ord->name);
    452536            psFree(values);
    453537            return false;
    454538        }
     539        if (!isfinite(value)) {
     540            psError(PS_ERR_UNKNOWN, true, "Value for DARK.ORDINATE %s is NAN", ord->name);
     541            psFree(values);
     542            return false;
     543        }
    455544        values->data.F32[i] = value;
    456545    }
    457     float norm = NAN;                   // Normalisation value
     546    double norm = NAN;                   // Normalisation value
    458547    bool doNorm = false;                // Do normalisation?
    459548    if (normConcept && strlen(normConcept) > 0) {
    460         if (!ordinateLookup(&norm, &inRange, normConcept, false, NAN, NAN, readout)) {
     549        if (!ordinateLookup(&norm, &inRange, normConcept, NULL, false, NAN, NAN, readout)) {
    461550            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept);
    462551            psFree(values);
    … …  
    470559        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information
    471560        if (ord->order <= 0) {
    472             psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order);
     561            psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order);
    473562            psFree(values);
    474563            psFree(orders);
    … …  
    701790            pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
    702791            psMetadata *row = psMetadataAlloc(); // FITS table row
    703             psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
     792            psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "DARK.ORDINATE name", ord->name);
     793
     794            // XXX write a dummy value if ord->rule == NULL? (eg, NONE)
     795            if (ord->rule) {
     796                psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", ord->rule);
     797            } else {
     798                psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", "NONE");
     799            }
     800
    704801            psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
    705802            psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
    … …  
    878975              ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);
    879976
     977              // load the ordinate rule; it is not an error if this field is missing or NULL
     978              // a NULL rule means 'use the name as the single concept'
     979              const char *rule = psMetadataLookupStr(&mdok, row, PM_DARK_FITS_RULE);
     980              if (!rule || !strcasecmp(rule, "NONE")) {
     981                  ord->rule = NULL;
     982              } else {
     983                  ord->rule = psStringCopy(rule);
     984              }
    880985              ordinates->data[i] = ord;
    881986          }
    … …  
    892997}
    893998
     999#if (HAVE_KAPA)
     1000#include <kapa.h>
     1001#include "pmKapaPlots.h"
     1002#include "pmVisual.h"
     1003
     1004static int nKapa = 0;
     1005static int *kapa = NULL;
     1006static bool plotFlag = true;
     1007static psArray *xVectors = NULL;
     1008
     1009// this init function only gets the ordinates for the first readout...
     1010bool pmDarkVisualInit(psArray *values) {
     1011   
     1012    if (!pmVisualIsVisual()) return true;
     1013
     1014    // skip if we have already opened the windows (or if none are requested...)
     1015    if (nKapa) return true;
     1016
     1017    // values has Ninput vectors with Norder elements; we need Norder vectors with Ninput elements...
     1018    int nOrders = 0;
     1019    for (int i = 0; i < values->n; i++) {
     1020        psVector *vect = values->data[i];
     1021        if (!nOrders) {
     1022            nOrders = vect->n;
     1023        } else {
     1024            psAssert (nOrders == vect->n, "mismatch in order vector lengths");
     1025        }
     1026    }
     1027    xVectors = psArrayAlloc(nOrders);
     1028    for (int i = 0; i < nOrders; i++) {
     1029        xVectors->data[i] = psVectorAlloc(values->n, PS_TYPE_F32);
     1030    }
     1031
     1032    for (int i = 0; i < values->n; i++) {
     1033        psVector *vect = values->data[i];
     1034        for (int j = 0; j < vect->n; j++) {
     1035            psVector *xVec = xVectors->data[j];
     1036            xVec->data.F32[i] = vect->data.F32[j];
     1037        }
     1038    }
     1039
     1040    nKapa = nOrders;
     1041
     1042    kapa = psAlloc(nKapa*sizeof(int));
     1043   
     1044    for (int i = 0; i < nKapa; i++) {
     1045        kapa[i] = -1;
     1046        pmVisualInitWindow(&kapa[i], "ppmerge");
     1047    }
     1048    return true;
     1049}
     1050
     1051bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) {
     1052
     1053    Graphdata graphdata;
     1054
     1055    if (!pmVisualIsVisual()) return true;
     1056
     1057    KapaInitGraph(&graphdata);
     1058
     1059    psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n);
     1060
     1061    psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32);
     1062    psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32);
     1063
     1064    for (int i = 0; i < xVectors->n; i++) {
     1065        psVector *x = xVectors->data[i];
     1066
     1067        // generate vectors of the unmasked values
     1068        int nSub = 0;
     1069        for (int j = 0; j < pixels->n; j++) {
     1070            if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue;
     1071            xSub->data.F32[nSub] = x->data.F32[j];
     1072            ySub->data.F32[nSub] = pixels->data.F32[j];
     1073            nSub ++;
     1074        }
     1075        xSub->n = ySub->n = nSub;
     1076       
     1077        // plot the unmasked values
     1078        pmVisualScaleGraphdata (&graphdata, xSub, ySub, false);
     1079        KapaSetGraphData(kapa[i], &graphdata);
     1080        KapaSetLimits(kapa[i], &graphdata);
     1081        KapaClearPlots (kapa[i]);
     1082
     1083        KapaSetFont (kapa[i], "courier", 14);
     1084        KapaBox (kapa[i], &graphdata);
     1085        KapaSendLabel (kapa[i], "ordinate", KAPA_LABEL_XM);
     1086        KapaSendLabel (kapa[i], "pixel values", KAPA_LABEL_YM);
     1087
     1088        graphdata.color = KapaColorByName("black");
     1089        graphdata.style = 2;
     1090        graphdata.ptype = 2;
     1091        KapaPrepPlot  (kapa[i], xSub->n, &graphdata);
     1092        KapaPlotVector(kapa[i], xSub->n, xSub->data.F32, "x");
     1093        KapaPlotVector(kapa[i], xSub->n, ySub->data.F32, "y");
     1094    }
     1095    pmVisualAskUser (&plotFlag);
     1096    return true;
     1097}
     1098
     1099bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) {
     1100
     1101    Graphdata graphdata;
     1102
     1103    if (!pmVisualIsVisual()) return true;
     1104
     1105    KapaInitGraph(&graphdata);
     1106
     1107    psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n);
     1108
     1109    psVector *yFit = psVectorAlloc(values->n, PS_TYPE_F32);
     1110
     1111    for (int i = 0; i < values->n; i++) {
     1112        psVector *coord = values->data[i];
     1113        yFit->data.F32[i] = psPolynomialMDEval (poly, coord);
     1114    }
     1115
     1116    for (int i = 0; i < xVectors->n; i++) {
     1117        psVector *xFit = xVectors->data[i];
     1118
     1119        KapaGetGraphData(kapa[i], &graphdata);
     1120        graphdata.color = KapaColorByName("red");
     1121        graphdata.style = 2;
     1122        graphdata.ptype = 7;
     1123        KapaPrepPlot  (kapa[i], xFit->n, &graphdata);
     1124        KapaPlotVector(kapa[i], xFit->n, xFit->data.F32, "x");
     1125        KapaPlotVector(kapa[i], xFit->n, yFit->data.F32, "y");
     1126    }
     1127    pmVisualAskUser (&plotFlag);
     1128    return true;
     1129}
     1130
     1131bool pmDarkVisualCleanup() {
     1132
     1133    for (int i = 0; i < nKapa; i++) {
     1134        KapaClose(kapa[i]);
     1135    }
     1136    psFree (kapa);
     1137    psFree (xVectors);
     1138    return true;
     1139}
     1140
     1141# else
     1142
     1143bool pmDarkVisualInit(psArray *values) { return true; }
     1144bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) { return true; }
     1145bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) { return true; }
     1146bool pmDarkVisualCleanup() { return true; }
     1147
     1148# endif

  • branches/pap/psModules/src/detrend/pmDark.h

    r21183 r25027  
    1313// An ordinate for fitting darks
    1414typedef struct {
    15     psString name;                      // Name of concept to fit
     15    psString name;                      // Name of ordinate
     16    psString rule;                      // Rule for generating ordinate (math on concepts)
    1617    int order;                          // Polynomial order to fit
    1718    bool scale;                         // Rescale values?
    … …  
    116117    );
    117118
     119bool pmDarkVisualInit(psArray *values);
     120bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask);
     121bool pmDarkVisualCleanup();
     122bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values);
    118123
    119124#endif

  • branches/pap/psModules/src/detrend/pmDetrendDB.c

    r23600 r25027  
    105105        DETREND_STRING_CASE(FRINGE);
    106106        DETREND_STRING_CASE(ASTROM);
     107        DETREND_STRING_CASE(NOISEMAP);
    107108    default:
    108109        return NULL;

  • branches/pap/psModules/src/detrend/pmDetrendDB.h

    r23487 r25027  
    3535    PM_DETREND_TYPE_BACKGROUND,
    3636    PM_DETREND_TYPE_ASTROM,
     37    PM_DETREND_TYPE_NOISEMAP,
    3738} pmDetrendType;
    3839

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

    r21457 r25027  
    6060    }
    6161
     62    // NOISEMAP : for now, not applied in the threaded loop
     63
    6264    return true;
    6365}

  • branches/pap/psModules/src/detrend/pmFringeStats.c

    r23259 r25027  
    464464    for (long i = 0; i < fringes->n; i++) {
    465465        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
     466        if (!fringe) {
     467            psWarning ("skipping empty fringe stats for entry %ld -- video cell?\n", i);
     468            continue;
     469        }
    466470        pmFringeRegions *regions = fringe->regions; // The fringe regions
    467471        if (numPoints == 0) {
    … …  
    496500    for (long i = 0; i < fringes->n; i++) {
    497501        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
     502        if (!fringe) {
     503            psWarning ("skipping empty fringe stats for entry %ld -- video cell?\n", i);
     504            continue;
     505        }
    498506        pmFringeRegions *regions = fringe->regions; // The fringe regions
    499507        // Copy the data over
    … …  
    520528//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    521529
    522 bool pmFringesFormat(pmCell *cell, psMetadata *header, const psArray *fringes)
    523 {
    524     PS_ASSERT_PTR_NON_NULL(cell, false);
     530psArray *pmFringesFormatTable(psMetadata *header, const psArray *fringes)
     531{
     532    PS_ASSERT_PTR_NON_NULL(header, false);
    525533    PS_ASSERT_ARRAY_NON_NULL(fringes, false);
    526534
    … …  
    531539        if (stats->regions != regions) {
    532540            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Regions for fringe statistics are not identical.\n");
    533             return false;
     541            return NULL;
    534542        }
    535543    }
    … …  
    557565    // Vectors: x, y, mask, f, df
    558566
    559     psMetadata *scalars = psMemIncrRefCounter(header); // Metadata to hold the scalars; will be the header
    560     if (!scalars) {
    561         scalars = psMetadataAlloc();
    562     }
    563     psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
    564                      regions->dX);
    565     psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height",
    566                      regions->dY);
    567     psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x",
    568                      regions->nX);
    569     psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y",
    570                      regions->nY);
    571     psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.HEADER", PS_DATA_METADATA,
    572                   "Header for fringe data", scalars);
    573     psFree(scalars);
    574 
     567    psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width", regions->dX);
     568    psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height", regions->dY);
     569    psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x", regions->nX);
     570    psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y", regions->nY);
    575571
    576572    psArray *table = psArrayAlloc(numRows); // The table
    … …  
    605601    }
    606602
    607     psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE", PS_DATA_ARRAY, "Fringe data", table);
    608 
    609     return true;
    610 }
    611 
    612 
    613 psArray *pmFringesParse(pmCell *cell)
    614 {
    615     PS_ASSERT_PTR_NON_NULL(cell, NULL);
     603    return table;
     604}
     605
     606psArray *pmFringesParseTable(psArray *table, psMetadata *header)
     607{
     608    PS_ASSERT_PTR_NON_NULL(table, NULL);
     609    PS_ASSERT_PTR_NON_NULL(header, NULL);
    616610
    617611    bool mdok;                          // Status of MD lookup
    618     psMetadata *header = psMetadataLookupMetadata(&mdok, cell->analysis, "FRINGE.HEADER"); // Header
    619     if (!mdok || !header) {
    620         psError(PS_ERR_UNKNOWN, false, "Unable to find header for fringe data.\n");
    621         return NULL;
    622     }
    623 
    624     psArray *table = psMetadataLookupPtr(NULL, cell->analysis, "FRINGE"); // FITS table
    625     if (!table) {
    626         psError(PS_ERR_UNKNOWN, false, "Unable to find table for fringe data.\n");
    627         return NULL;
    628     }
    629 
    630612
    631613    // Read the scalars from the header
    … …  
    722704}
    723705
     706bool pmFringesFormat(pmCell *cell, psMetadata *inHeader, const psArray *fringes)
     707{
     708    PS_ASSERT_PTR_NON_NULL(cell, false);
     709    PS_ASSERT_ARRAY_NON_NULL(fringes, false);
     710
     711    psMetadata *header = psMemIncrRefCounter(inHeader); // Metadata to hold the scalars; will be the header
     712    if (!inHeader) {
     713        inHeader = psMetadataAlloc();
     714    }
     715
     716    psArray *table = pmFringesFormatTable(header, fringes);
     717    if (!table) {
     718        psError (PS_ERR_UNKNOWN, false, "unable to generate table from fringes");
     719        psFree(header);
     720        return NULL;
     721    }
     722
     723    psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.HEADER", PS_DATA_METADATA, "Header for fringe data", header);
     724    psFree(header);
     725
     726    psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.TABLE", PS_DATA_ARRAY, "Fringe data", table);
     727    psFree(table);
     728
     729    return true;
     730}
     731
     732psArray *pmFringesParse(pmCell *cell)
     733{
     734    PS_ASSERT_PTR_NON_NULL(cell, NULL);
     735
     736    bool mdok;                          // Status of MD lookup
     737    psMetadata *header = psMetadataLookupMetadata(&mdok, cell->analysis, "FRINGE.HEADER"); // Header
     738    if (!mdok || !header) {
     739        psError(PS_ERR_UNKNOWN, false, "Unable to find header for fringe data.\n");
     740        return NULL;
     741    }
     742
     743    psArray *table = psMetadataLookupPtr(NULL, cell->analysis, "FRINGE.TABLE"); // FITS table
     744    if (!table) {
     745        psError(PS_ERR_UNKNOWN, false, "Unable to find table for fringe data.\n");
     746        return NULL;
     747    }
     748
     749    psArray *fringes = pmFringesParseTable(table, header);
     750    if (!fringes) {
     751        psError(PS_ERR_UNKNOWN, false, "Unable to add parse the fringe table data");
     752        return NULL;
     753    }
     754
     755    return fringes;
     756}
    724757
    725758//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    821854
    822855    // Solve the least-squares equation
    823     if (! psMatrixGJSolve(A, B)) {
    824         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
     856    if (!psMatrixGJSolve(A, B)) {
     857        psLogMsg("psModules.detrend", PS_LOG_INFO, "Could not solve linear equations.  Returning NULL.\n");
     858        psFree(A);
     859        psFree(B);
    825860        return false;
    826861    }
    … …  
    882917
    883918    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics
    884     psVectorStats(stats, diffs, NULL, mask, 1);
     919    if (!psVectorStats(stats, diffs, NULL, mask, 1)) {
     920        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     921        return 0;
     922    }
    885923    float middle = stats->sampleMedian; // The middle of the distribution
    886924    float thresh = rej * 0.74 * (stats->sampleUQ - stats->sampleLQ); // The rejection threshold
    … …  
    9691007
    9701008    // Get rid of the extreme outliers by assuming most of the points are somewhat clustered
    971     psVectorStats(median, science->f, NULL, NULL, 0);
     1009    if (!psVectorStats(median, science->f, NULL, NULL, 0)) {
     1010        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1011        return NULL;
     1012    }
    9721013    scale->coeff->data.F32[0] = median->sampleMedian;
    9731014    for (int i = 0; i < fringes->n; i++) {
    9741015        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
    975         psVectorStats(median, fringe->f, NULL, NULL, 0);
     1016        if (!psVectorStats(median, fringe->f, NULL, NULL, 0)) {
     1017            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1018            return NULL;
     1019        }
    9761020        scale->coeff->data.F32[0] -= median->sampleMedian;
    9771021        scale->coeff->data.F32[i] = 0.0;

  • branches/pap/psModules/src/detrend/pmFringeStats.h

    r21183 r25027  
    142142//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    143143
    144 /// Write an array of fringes measurements to a FITS table.
    145 ///
    146 /// Writes an array of fringe measurements for a cell as a FITS table in the analysis metadata.  The array of
    147 /// fringe statistics must all use the same fringe regions (or there is no point in storing them all
    148 /// together).  The header is supplemented with scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY,
    149 /// PSFRNGNX, PSFRNGNY) from the fringe regions, while the fringe coordinates and mask are written as a FITS
    150 /// table (as x, y, mask, f, df; f and df are vectors).
     144/// Convert an array of fringes measurements to a psArray suitable for writing as a FITS table.
     145///
     146/// Converts an array of fringe measurements for a cell into the corresponding rows of a FITS
     147/// table (array of psMetadata).  The array of fringe statistics must all use the same fringe
     148/// regions (or there is no point in storing them all together).  The header is supplemented
     149/// with scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY, PSFRNGNX, PSFRNGNY) from the
     150/// fringe regions, while the fringe coordinates and mask are written as a FITS table rows (as
     151/// x, y, mask, f, df; f and df are vectors).  Use psFitsTableWrite to save the resulting rows
     152/// to disk.
     153psArray *pmFringesFormatTable(psMetadata *header, const psArray *fringes);
     154
     155
     156/// Parses an array of fringes measurements from a FITS table.
     157///
     158/// The fringes for the cell are read from the FITS table (array of psMetadata rows).  The
     159/// table provides the region and the (possibly multiple) fringe statistics for that region.
     160/// The supplied header defines the scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY,
     161/// PSFRNGNX, PSFRNGNY)
     162psArray *pmFringesParseTable(psArray *table, psMetadata *header);
     163
     164/// Deprecated Function: converts the fringes to a FITS table (array of psMetadata) and saves
     165/// them on the cell->analysis; scalar values are dX, dY, nX and nY are written to the header
     166/// (as PSFRNGDX, PSFRNGDY, PSFRNGNX, PSFRNGNY)
    151167bool pmFringesFormat(pmCell *cell,   ///< Cell for which to write
    152168                     psMetadata *header, ///< Header, or NULL
    … …  
    154170                    );
    155171
    156 /// Parses an array of fringes measurements from a FITS table.
    157 ///
    158 /// The fringes for the cell are read from the FITS table in the analysis metadata.  The table provides the
    159 /// region and the (possibly multiple) fringe statistics for that region.  The current extension is used if
    160 /// the extension name is not provided.
     172/// Deprecated Function: pulls a the header and FITS table (array of psMetadata) representing
     173/// the fringes from the cell->analysis and converts to fringe measurements
    161174psArray *pmFringesParse(pmCell *cell ///< Cell for which to read fringes
    162175                       );
    163 
    164176
    165177//////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • branches/pap/psModules/src/detrend/pmMaskBadPixels.c

    r21183 r25027  
    8181
    8282
    83 bool pmMaskFlagSuspectPixels(pmReadout *output, const pmReadout *readout, float median, float stdev,
    84                              float rej, psImageMaskType maskVal)
     83bool pmMaskFlagSuspectPixelsBySigma(pmReadout *output, const pmReadout *readout, float median, float stdev,
     84                                    float rej, psImageMaskType maskVal)
    8585{
    8686    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    115115        // If we get down here and the statistics are missing, then we should go and mask the entire image
    116116        psWarning("Missing statistics --- flagging entire image as suspect.");
    117         return (psImage*)psBinaryOp(suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_F32));
     117        psBinaryOp (suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_F32));
     118        return true;
    118119    }
    119120
    … …  
    128129        for (int x = 0; x < image->numCols; x++) {
    129130            if (fabs((image->data.F32[y][x] - median) / stdev) < rej) continue;
     131            if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) continue;
     132            suspect->data.F32[y][x] += 1.0;
     133        }
     134    }
     135    psFree(suspect);                    // Drop reference
     136
     137    psMetadataItem *numItem = psMetadataLookup(output->analysis, PM_MASK_ANALYSIS_NUM); // Item with number
     138    assert(numItem);
     139    numItem->data.S32++;
     140
     141    return true;
     142}
     143
     144bool pmMaskFlagSuspectPixelsByValue(pmReadout *output, const pmReadout *readout,
     145                                    float min, float max, psImageMaskType maskVal)
     146{
     147    PS_ASSERT_PTR_NON_NULL(readout, false);
     148    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
     149    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false);
     150    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
     151    if (readout->mask) {
     152        PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false);
     153        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false);
     154        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_IMAGE_MASK, false);
     155    }
     156    PS_ASSERT_PTR_NON_NULL(output, false);
     157
     158    bool mdok;                          // Status of MD lookup
     159    psImage *suspect = psMetadataLookupPtr(&mdok, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img
     160    if (suspect) {
     161        PS_ASSERT_IMAGE_NON_EMPTY(suspect, false);
     162        PS_ASSERT_IMAGE_TYPE(suspect, PS_TYPE_F32, false);
     163        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, suspect, false);
     164        psMemIncrRefCounter(suspect);
     165    } else {
     166        suspect = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
     167        psImageInit(suspect, 0);
     168        psMetadataAddImage(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_SUSPECT, PS_META_REPLACE,
     169                           "Suspect pixels", suspect);
     170        psMetadataAddS32(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_NUM, PS_META_REPLACE,
     171                         "Number of input images", 0);
     172    }
     173
     174    psImage *image = readout->image;    // Image of interest
     175    psImage *mask = readout->mask;      // Corresponding mask
     176
     177    psTrace ("psModules.detrend", 3, "suspect: < %f or > %f\n", min, max);
     178
     179    // XXX this loop could result in pixels with suspect = 0.0 but no valid input pixels (all
     180    // masked).  need to track the number of good as well as suspect pixels?
     181    for (int y = 0; y < image->numRows; y++) {
     182        for (int x = 0; x < image->numCols; x++) {
     183            bool above = image->data.F32[y][x] > max;
     184            bool below = image->data.F32[y][x] < min;
     185            if (!above && !below) continue;
    130186            if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) continue;
    131187            suspect->data.F32[y][x] += 1.0;

  • branches/pap/psModules/src/detrend/pmMaskBadPixels.h

    r21183 r25027  
    5050/// high value in the suspect pixels image, allowing them to be identified.  The suspect pixels
    5151/// image is of type S32.  The relevant median and standard deviation must be supplied in the
    52 /// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEVe
    53 bool pmMaskFlagSuspectPixels(pmReadout *output, ///< Output readout, optionally with suspect pixels image
     52/// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEV
     53bool pmMaskFlagSuspectPixelsBySigma(pmReadout *output, ///< Output readout, optionally with suspect pixels image
    5454                             const pmReadout *readout, ///< Readout to inspect
    5555                             float median, ///< Image median
    5656                             float stdev, ///< Image standard deviation
    5757                             float rej, ///< Rejection threshold (standard deviations)
     58                             psImageMaskType maskVal ///< Mask value for statistics
     59    );
     60
     61/// Find out-of-range pixels and flag them
     62///
     63/// Pixels great > max or < min have the corresponding pixel in the "suspect pixels" image
     64/// incremented.  After accumulating over a suitable sample of images, bad pixels should have a
     65/// high value in the suspect pixels image, allowing them to be identified.  The suspect pixels
     66/// image is of type S32.  The relevant median and standard deviation must be supplied in the
     67/// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEV
     68bool pmMaskFlagSuspectPixelsByValue(pmReadout *output, ///< Output readout, optionally with suspect pixels image
     69                             const pmReadout *readout, ///< Readout to inspect
     70                             float min, ///< Image min acceptable value
     71                             float max, ///< Image max acceptable value
    5872                             psImageMaskType maskVal ///< Mask value for statistics
    5973    );

  • branches/pap/psModules/src/detrend/pmOverscan.c

    r23826 r25027  
    7474            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    7575            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
    76             psVectorStats(myStats, values, NULL, NULL, 0);
     76            if (!psVectorStats(myStats, values, NULL, NULL, 0)) {
     77                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     78                return false;
     79            }
    7780            reduced->data.F32[i] = psStatsGetValue(myStats, statistic);
    7881        } else if (overscanOpts->fitType == PM_FIT_NONE) {
    … …  
    275278        psFree(iter);
    276279
    277         (void)psVectorStats(stats, pixels, NULL, NULL, 0);
     280        if (!psVectorStats(stats, pixels, NULL, NULL, 0)) {
     281            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     282            return false;
     283        }
    278284        psFree(pixels);
    279285        double reduced = psStatsGetValue(stats, statistic); // Result of statistics
    … …  
    349355          psString comment = NULL;    // Comment to add
    350356          psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    351           psVectorStats (vectorStats, reduced, NULL, NULL, 0);
     357          if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
     358              psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     359              return false;
     360          }
    352361          psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
    353362          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    … …  
    412421          psString comment = NULL;    // Comment to add
    413422          psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    414           psVectorStats (vectorStats, reduced, NULL, NULL, 0);
     423          if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
     424              psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     425              return false;
     426          }
    415427          psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
    416428          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");

  • branches/pap/psModules/src/detrend/pmRemnance.c

    r23259 r25027  
    6666            values->n = numValues;
    6767            if (!psVectorStats(stats, values, NULL, NULL, 0)) {
    68                 // Can't do anything about it
    69                 psErrorClear();
     68                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     69                return false;
     70            }
     71            if (isnan(stats->sampleMedian)) {
    7072                maxMask = max;
    7173                continue;

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

    r23487 r25027  
    350350    psStats *rawStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    351351    psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    352     psVectorStats (rawStats, counts, NULL, NULL, 0);
    353     psVectorStats (resStats, resid, NULL, NULL, 0);
     352    if (!psVectorStats (rawStats, counts, NULL, NULL, 0)) {
     353        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     354        return NULL;
     355    }
     356    if (!psVectorStats (resStats, resid, NULL, NULL, 0)) {
     357        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     358        return NULL;
     359    }
    354360
    355361    // XXX temporary hard-wired minimum stdev improvement factor
    356362    psTrace("psModules.detrend", 3, "raw scatter %f vs res scatter %f\n", rawStats->sampleStdev, resStats->sampleStdev);
    357363    if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false;
     364    if (isnan(rawStats->sampleStdev) || isnan(resStats->sampleStdev)) corr->valid = false;
    358365
    359366    psFree (rawStats);

  • branches/pap/psModules/src/detrend/pmSkySubtract.c

    r21183 r25027  
    410410    // XXX: How do we know if these matrix operations were successful?
    411411    //
    412     Aout = psMatrixLUD(Aout, &outPerm, A);
     412    Aout = psMatrixLUDecomposition(Aout, &outPerm, A);
    413413    PS_ASSERT_IMAGE_NON_NULL(Aout, NULL);
    414414    PS_ASSERT_IMAGE_NON_EMPTY(Aout, NULL);
    415415    psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
    416     psMatrixLUSolve(C, Aout, B, outPerm);
     416    psMatrixLUSolution(C, Aout, B, outPerm);
    417417
    418418    //

  • branches/pap/psModules/src/extras/Makefile.am

    r21422 r25027  
    88        psIOBuffer.c \
    99        pmKapaPlots.c \
    10         pmVisual.c
     10        pmVisual.c \
     11        ippStages.c
    1112
    1213pkginclude_HEADERS = \
    … …  
    1516        psIOBuffer.h \
    1617        pmKapaPlots.h \
    17         pmVisual.h
     18        pmVisual.h \
     19        ippStages.h
    1820
    1921CLEANFILES = *~

  • branches/pap/psModules/src/extras/pmVisual.c

    r23242 r25027  
    281281    psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    282282    psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    283     psVectorStats (statsX, xVec, NULL, NULL, 0);
    284     psVectorStats (statsY, yVec, NULL, NULL, 0);
     283    if (!psVectorStats (statsX, xVec, NULL, NULL, 0)) {
     284        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     285        return false;
     286    }
     287    if (!psVectorStats (statsY, yVec, NULL, NULL, 0)) {
     288        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     289        return false;
     290    }
    285291
    286292    float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;

  • branches/pap/psModules/src/imcombine/pmImageCombine.c

    r21183 r25027  
    132132        // Combine all the pixels, using the specified stat.
    133133        if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, 0xff)) {
     134            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     135            return false;
     136        }
     137        if (isnan(stats->sampleMean)) {
    134138            combine->data.F32[y][x] = NAN;
    135139            psFree(buffer);
    … …  
    137141        }
    138142        float combinedPixel = stats->sampleMean; // Value of the combination
    139 
     143       
    140144        if (iter == 0) {
    141145            // Save the value produced with no rejection, since it may be useful later
    … …  
    364368    // Get the median
    365369    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    366     psVectorStats(stats, pixels, NULL, mask, 1);
     370    if (!psVectorStats(stats, pixels, NULL, mask, 1)) {
     371        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     372    }
    367373    float median = stats->sampleMedian;
    368374    psFree(stats);

  • branches/pap/psModules/src/imcombine/pmPSFEnvelope.c

    r23242 r25027  
    3333
    3434
    35 //#define TESTING                         // Enable test output
     35// #define TESTING                         // Enable test output
    3636#define PEAK_FLUX 1.0e4                 // Peak flux for each source
    3737#define SKY_VALUE 0.0e0                 // Sky value for fake image
    … …  
    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
     42#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
     43                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
    4244
    4345
    … …  
    112114    psImageInit(envelope, SKY_VALUE);
    113115    pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
    114     float maxRadius = 0.0;              // Maximum radius of sources
     116    float maxRadius = 0.0;              // Maximum radius for sources
     117    psVector *numbers = psVectorAlloc(numFakes, PS_TYPE_S32); // Number of detections for each source
     118    psVectorInit(numbers, 0);
    115119    for (int i = 0; i < inputs->n; i++) {
    116120        pmPSF *psf = inputs->data[i];   // PSF of interest
    … …  
    127131            psFree(xOffset);
    128132            psFree(fakes);
     133            psFree(numbers);
    129134            psf->residuals = resid;
    130135            return NULL;
    … …  
    140145
    141146            double flux = fakeRO->image->data.F32[(int)y][(int)x];
     147            if (!isfinite(flux) || flux < 0) {
     148                continue;
     149            }
    142150            float norm = PEAK_FLUX / flux; // Normalisation for source
    143151            psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region
    … …  
    151159            // Get the radius
    152160            pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source
    153             psAssert (model, "failed to generate model: should this be an error or not?");
     161            if (!model || (model->flags & MODEL_MASK)) {
     162                continue;
     163            }
    154164            float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source
     165            if (srcRadius == 0) {
     166                continue;
     167            }
    155168            if (srcRadius > maxRadius) {
    156169                maxRadius = srcRadius;
    157170            }
    158             psFree(model);
     171
     172            // If we got this far, the source is decent
     173            numbers->data.S32[j]++;
    159174        }
    160175
    … …  
    175190    psFree(fakeRO);
    176191
    177     if (maxRadius > radius) {
    178         maxRadius = radius;
    179     }
    180 
    181192#ifdef TESTING
    182193    {
    … …  
    190201
    191202    // Put the fake sources onto a full-size image
     203    psArray *goodFakes = psArrayAllocEmpty(numFakes); // Good fake sources
    192204    pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels
    193205    readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    … …  
    195207    for (int i = 0; i < numFakes; i++) {
    196208        pmSource *source = fakes->data[i]; // Fake source
     209        if (numbers->data.S32[i] > 0) {
     210            psArrayAdd(goodFakes, goodFakes->n, source);
     211        }
     212
    197213        // Position of source on fake image
    198214        int xFake = source->peak->x + xOffset->data.S32[i];
    … …  
    213229            psFree(yOffset);
    214230            psFree(fakes);
     231            psFree(numbers);
    215232            return NULL;
    216233        }
    … …  
    220237    psFree(yOffset);
    221238    psFree(envelope);
     239    psFree(numbers);
     240
     241    psFree(fakes);
     242    fakes = goodFakes;
     243    numFakes = fakes->n;
     244
     245    if (numFakes == 0.0) {
     246        psError(PS_ERR_UNKNOWN, false, "No fake sources are suitable for PSF fitting.");
     247        psFree(fakes);
     248        psFree(readout);
     249        return false;
     250    }
    222251
    223252    // XXX Setting the variance seems to be an art
    … …  
    232261    psImageInit(readout->mask, 0);
    233262
     263    if (maxRadius > radius) {
     264        maxRadius = radius;
     265    }
     266
    234267#ifdef TESTING
    235268    {
    … …  
    243276
    244277    // Reset the sources to point to the new pixels, and measure the moments in preparation for PSF fitting
     278    int numMoments = 0;                 // Number of moments measured
    245279    for (int i = 0; i < numFakes; i++) {
    246280        pmSource *source = fakes->data[i]; // Fake source
    … …  
    257291        source->maskObj = NULL;
    258292
    259         if (!pmSourceDefinePixels(source, readout, x, y, maxRadius)) {
     293        if (!pmSourceDefinePixels(source, readout, x, y, radius)) {
    260294            psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
    261295            psFree(readout);
    … …  
    264298        }
    265299
    266         if (!pmSourceMoments(source, maxRadius)) {
    267             psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source.");
    268             psFree(readout);
    269             psFree(fakes);
    270             return NULL;
    271         }
     300        // measure the source moments: tophat windowing, no pixel S/N cutoff
     301        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
     302            // Can't do anything about it; limp along as best we can
     303            psErrorClear();
     304            continue;
     305        }
     306        numMoments++;
     307    }
     308
     309    if (numMoments == 0) {
     310        psError(PS_ERR_UNKNOWN, true, "Unable to measure moments for sources.");
     311        psFree(fakes);
     312        psFree(readout);
     313        return NULL;
    272314    }
    273315
    … …  
    286328    options->psfFieldXo = 0;
    287329    options->psfFieldYo = 0;
     330    options->chiFluxTrend = false;      // All sources have similar flux, so fitting a trend often fails
    288331
    289332    pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);

  • branches/pap/psModules/src/imcombine/pmReadoutCombine.c

    r21363 r25027  
    173173    for (int i = 0; i < inputs->n; i++) {
    174174        pmReadout *readout = inputs->data[i]; // Readout of interest
     175        psAssert(readout, "readout was not defined");
     176        if (!readout->process) continue;
    175177        if (!readout->image) {
    176178            psError(PS_ERR_UNEXPECTED_NULL, true, "Image data is missing for image %d.\n", i);
    … …  
    271273        for (int r = 0; r < inputs->n; r++) {
    272274            pmReadout *readout = inputs->data[r]; // Input readout
     275            if (!readout->process) continue;
    273276            psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r,
    274277                    minInputCols - readout->col0, maxInputCols - readout->col0,
    … …  
    277280    }
    278281    #endif
     282
     283    // set up windows for visualization (if selected)
     284    pmReadoutCombineVisualInit();
    279285
    280286    // Dereference output products
    … …  
    308314            for (int r = 0; r < inputs->n; r++) {
    309315                pmReadout *readout = inputs->data[r]; // Input readout
     316                if (!readout->process) {
     317                    maskData[r] = 1;
     318                    continue;
     319                }
    310320                int yIn = i - readout->row0; // y position on input readout
    311321                int xIn = j - readout->col0; // x position on input readout
    … …  
    384394            // Combination
    385395            if (!psVectorStats(stats, pixels, errors, mask, 1)) {
    386                 // Can't do much about it, but it's not worth worrying about
    387                 psErrorClear();
     396                psError(PS_ERR_UNKNOWN, false, "error in pixel stats");
     397                return false;
     398            }
     399           
     400            outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
     401
     402            if (!isfinite(outputImage[yOut][xOut])) {
     403                pmReadoutCombineVisualPixels(pixels, mask, outputImage[yOut][xOut]);
     404            }
     405
     406            if (isnan(outputImage[yOut][xOut])) {
    388407                outputImage[yOut][xOut] = NAN;
    389408                outputMask[yOut][xOut] = params->blank;
     409                sigma->data.F32[yOut][xOut] = NAN;
    390410                if (params->variances) {
    391411                    outputVariance[yOut][xOut] = NAN;
    392412                }
    393                 sigma->data.F32[yOut][xOut] = NAN;
    394             } else {
    395                 outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
    396                 outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank;
    397                 if (params->variances) {
    398                     float stdev = psStatsGetValue(stats, combineStdev);
    399                     outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance
    400                     // XXXX this is not the correct formal error.
    401                     // also, the weighted mean is not obviously the correct thing here
    402                 }
    403                 sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev);
    404             }
     413                continue;
     414            }
     415            outputMask[yOut][xOut] = 0;
     416            sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev);
     417            if (params->variances) {
     418                float stdev = psStatsGetValue(stats, combineStdev);
     419                outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance
     420                // XXXX this is not the correct formal error.
     421                // also, the weighted mean is not obviously the correct thing here
     422            }
    405423        }
    406424    }
    … …  
    423441}
    424442
     443#if (HAVE_KAPA)
     444#include <kapa.h>
     445#include "pmKapaPlots.h"
     446#include "pmVisual.h"
     447
     448static int kapa = -1;
     449static bool plotFlag = true;
     450
     451// this init function only gets the ordinates for the first readout...
     452bool pmReadoutCombineVisualInit(void) {
     453   
     454    if (!pmVisualIsVisual()) return true;
     455
     456    // skip if we have already opened the windows (or if none are requested...)
     457    if (kapa != -1) return true;
     458
     459    pmVisualInitWindow(&kapa, "ppmerge");
     460    return true;
     461}
     462
     463bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) {
     464
     465    Graphdata graphdata;
     466    float xline[2], yline[2];
     467   
     468    if (!pmVisualIsVisual()) return true;
     469
     470    if (!plotFlag) return true;
     471
     472    KapaInitGraph(&graphdata);
     473
     474    psVector *xAll = psVectorAlloc(pixels->n, PS_TYPE_F32);
     475    psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32);
     476    psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32);
     477
     478    // generate vectors of the unmasked values
     479    int nSub = 0;
     480    for (int j = 0; j < pixels->n; j++) {
     481        xAll->data.F32[j] = j;
     482        if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue;
     483        xSub->data.F32[nSub] = j;
     484        ySub->data.F32[nSub] = pixels->data.F32[j];
     485        nSub ++;
     486    }
     487    xSub->n = ySub->n = nSub;
     488    xAll->n = pixels->n;
     489       
     490    xline[0] = 0;
     491    xline[1] = pixels->n;
     492    yline[0] = mean;
     493    yline[1] = mean;
     494
     495    // plot the unmasked values
     496    pmVisualScaleGraphdata (&graphdata, xAll, pixels, false);
     497    KapaSetGraphData(kapa, &graphdata);
     498    KapaSetLimits(kapa, &graphdata);
     499    KapaClearPlots (kapa);
     500
     501    KapaSetFont (kapa, "courier", 14);
     502    KapaBox (kapa, &graphdata);
     503    KapaSendLabel (kapa, "ordinate", KAPA_LABEL_XM);
     504    KapaSendLabel (kapa, "pixel values", KAPA_LABEL_YM);
     505
     506    graphdata.color = KapaColorByName("black");
     507    graphdata.style = 2;
     508    graphdata.ptype = 2;
     509    KapaPrepPlot  (kapa, xSub->n, &graphdata);
     510    KapaPlotVector(kapa, xSub->n, xSub->data.F32, "x");
     511    KapaPlotVector(kapa, xSub->n, ySub->data.F32, "y");
     512
     513    graphdata.color = KapaColorByName("red");
     514    graphdata.style = 2;
     515    graphdata.ptype = 7;
     516    KapaPrepPlot  (kapa, xAll->n, &graphdata);
     517    KapaPlotVector(kapa, xAll->n, xAll->data.F32, "x");
     518    KapaPlotVector(kapa, xAll->n, pixels->data.F32, "y");
     519
     520    graphdata.color = KapaColorByName("blue");
     521    graphdata.style = 0;
     522    graphdata.ptype = 7;
     523    KapaPrepPlot  (kapa, 2, &graphdata);
     524    KapaPlotVector(kapa, 2, xline, "x");
     525    KapaPlotVector(kapa, 2, yline, "y");
     526
     527    pmVisualAskUser (&plotFlag);
     528    return true;
     529}
     530
     531bool pmReadoutCombineVisualCleanup(void) {
     532
     533    if (!pmVisualIsVisual()) return true;
     534    if (kapa == -1) return true;
     535
     536    KapaClose(kapa);
     537    return true;
     538}
     539
     540# else
     541
     542bool pmReadoutCombineVisualInit(void) { return true; }
     543bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) { return true; }
     544bool pmReadoutCombineVisualCleanup(void) { return true; }
     545
     546# endif

  • branches/pap/psModules/src/imcombine/pmReadoutCombine.h

    r21363 r25027  
    4747                     );
    4848
     49bool pmReadoutCombineVisualInit(void);
     50bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean);
     51bool pmReadoutCombineVisualCleanup(void);
     52
    4953/// @}
    5054#endif

  • branches/pap/psModules/src/imcombine/pmSubtraction.c

    r23839 r25027  
    196196{
    197197    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
    198     int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    199198    for (int y = -footprint; y <= footprint; y++) {
    200         // Convolution with a delta function is just the value specified by the offset
    201         memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes);
     199        for (int x = -footprint; x <= footprint; x++) {
     200            convolved->kernel[y][x] = image->kernel[y - v][x - u];
     201        }
    202202    }
    203203    return convolved;
    204204}
    205205
    206 // Take a subset of an image
    207 static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved
    208                                            psRegion region, // Region of interest (with border)
    209                                            bool threaded // Are we running threaded?
    210                                            )
    211 {
    212     if (!image) {
    213         return NULL;
    214     }
    215     // XXX if (threaded) {
    216     // XXX     psMutexLock(image);
    217     // XXX }
    218     psImage *subset = psImageSubset(image, region); // Subset image, to return
    219     // XXX if (threaded) {
    220     // XXX     psMutexUnlock(image);
    221     // XXX }
    222     return subset;
    223 }
    224 
    225 // Free the subset of an image
    226 static inline void convolveSubsetFree(psImage *parent, // Parent image
    227                                       psImage *child, // Child (subset) image
    228                                       bool threaded // Are we running threaded?
    229                                       )
    230 {
    231     if (!child || !parent) {
    232         return;
    233     }
    234     // XXX if (threaded) {
    235     // XXX     psMutexLock(parent);
    236     // XXX }
    237     psFree(child);
    238     // XXX if (threaded) {
    239     // XXX     psMutexUnlock(parent);
    240     // XXX }
    241     return;
    242 }
    243206
    244207// Convolve an image using FFT
    … …  
    256219                                  region.y0 - size, region.y1 + size); // Add a border
    257220
    258     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    259 
    260     psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve
    261     psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask
     221    psImage *subImage = image ? psImageSubset(image, border) : NULL; // Subimage to convolve
     222    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask
    262223
    263224    psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
    264225
    265     convolveSubsetFree(image, subImage, threaded);
    266     convolveSubsetFree(mask, subMask, threaded);
     226    psFree(subImage);
     227    psFree(subMask);
    267228
    268229    // Now, we have to stick it in where it belongs
    … …  
    300261                                  region.y0 - size, region.y1 + size); // Add a border
    301262
    302     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    303 
    304     psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map
    305     psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
    306     psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
     263    psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map
     264    psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image
     265    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask
    307266
    308267    // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
    … …  
    310269    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
    311270
    312     convolveSubsetFree(variance, subVariance, threaded);
    313     convolveSubsetFree(sys, subSys, threaded);
    314     convolveSubsetFree(mask, subMask, threaded);
     271    psFree(subVariance);
     272    psFree(subSys);
     273    psFree(subMask);
    315274
    316275    // Now, we have to stick it in where it belongs
    … …  
    420379        if (box > 0) {
    421380            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
    422             bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    423 
    424381            psRegion region = psRegionSet(colMin - box, colMax + box,
    425382                                          rowMin - box, rowMax + box); // Region to convolve
    426383
    427             psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve
     384            psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve
    428385
    429386            psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad,
    430387                                                     -box, box, -box, box); // Convolved subtraction mask
    431388
    432             convolveSubsetFree(subMask, image, threaded);
     389            psFree(image);
    433390
    434391            psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
    … …  
    653610                  float value = 0.0;    // Value of convolved pixel
    654611                  int uMin = x - size, uMax = x + size; // Range for u
    655                   psF32 *xKernelData = xKernel->data.F32; // Kernel values
     612                  psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
    656613                  psF32 *imageData = &image->kernel[y][uMin]; // Image values
    657                   for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) {
     614                  for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
    658615                      value += *xKernelData * *imageData;
    659616                  }
    … …  
    668625                  float value = 0.0;    // Value of convolved pixel
    669626                  int vMin = y - size, vMax = y + size; // Range for v
    670                   psF32 *yKernelData = yKernel->data.F32; // Kernel values
     627                  psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
    671628                  psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
    672                   for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) {
     629                  for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
    673630                      value += *yKernelData * *imageData;
    674631                  }
    … …  
    821778    }
    822779    psFree(mask);
     780
     781    // XXX raise an error?
     782    if (isnan(stats->sampleMean)) {
     783        psFree(stats);
     784        return -1;
     785    }
    823786
    824787    double mean, rms;                 // Mean and RMS of deviations

  • branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

    r23839 r25027  
    1515#include "pmSubtractionVisual.h"
    1616
    17 //#define TESTING
     17// #define TESTING                         // TESTING output for debugging; may not work with threads!
     18
     19#define USE_VARIANCE                    // Include variance in equation?
    1820
    1921//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    3133    for (int y = - footprint; y <= footprint; y++) {
    3234        for (int x = - footprint; x <= footprint; x++) {
    33             sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x];
     35            double value = image1->kernel[y][x] * image2->kernel[y][x];
     36#ifdef USE_VARIANCE
     37            value /= variance->kernel[y][x];
     38#endif
     39            sum += value;
    3440        }
    3541    }
    … …  
    195201            for (int y = - footprint; y <= footprint; y++) {
    196202                for (int x = - footprint; x <= footprint; x++) {
    197                     sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     203                    double value = conv->kernel[y][x];
     204#ifdef USE_VARIANCE
     205                    value /= variance->kernel[y][x];
     206#endif
     207                    sumC += value;
    198208                }
    199209            }
    … …  
    218228        for (int y = - footprint; y <= footprint; y++) {
    219229            for (int x = - footprint; x <= footprint; x++) {
    220                 double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x];
     230                double invNoise2 = 1.0;
     231#ifdef USE_VARIANCE
     232                invNoise2 /= variance->kernel[y][x];
     233#endif
    221234                double value = input->kernel[y][x] * invNoise2;
    222235                sumI += value;
    … …  
    277290        for (int y = - footprint; y <= footprint; y++) {
    278291            for (int x = - footprint; x <= footprint; x++) {
    279                 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     292                double value = target->kernel[y][x] * conv->kernel[y][x];
     293#ifdef USE_VARIANCE
     294                value /= variance->kernel[y][x];
     295#endif
     296                sumTC += value;
    280297            }
    281298        }
    … …  
    297314        for (int y = - footprint; y <= footprint; y++) {
    298315            for (int x = - footprint; x <= footprint; x++) {
    299                 float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x];
     316                double value = target->kernel[y][x];
     317#ifdef USE_VARIANCE
     318                value /= variance->kernel[y][x];
     319#endif
    300320                sumIT += value * input->kernel[y][x];
    301321                sumT += value;
    … …  
    366386        for (int y = - footprint; y <= footprint; y++) {
    367387            for (int x = - footprint; x <= footprint; x++) {
    368                 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     388                double value = conv->kernel[y][x];
     389#ifdef USE_VARIANCE
     390                value /= variance->kernel[y][x];
     391#endif
     392                sumC += value;
    369393            }
    370394        }
    … …  
    765789        for (int i = 0; i < stamps->num; i++) {
    766790            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     791
    767792            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     793
     794#ifdef TESTING
     795              // XXX double-check for NAN in data:
     796                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
     797                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
     798                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
     799                            fprintf (stderr, "WARNING: NAN in matrix1\n");
     800                        }
     801                    }
     802                }
     803                for (int ix = 0; ix < stamp->vector1->n; ix++) {
     804                    if (!isfinite(stamp->vector1->data.F64[ix])) {
     805                        fprintf (stderr, "WARNING: NAN in vector1\n");
     806                    }
     807                }
     808#endif
     809
    768810                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
    769811                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
    … …  
    774816            }
    775817        }
     818
     819#ifdef TESTING
     820        for (int ix = 0; ix < sumVector->n; ix++) {
     821            if (!isfinite(sumVector->data.F64[ix])) {
     822                fprintf (stderr, "WARNING: NAN in vector1\n");
     823            }
     824        }
     825#endif
     826
    776827        calculatePenalty(sumVector, kernels);
    777828
    778829#ifdef TESTING
     830        for (int ix = 0; ix < sumVector->n; ix++) {
     831            if (!isfinite(sumVector->data.F64[ix])) {
     832                fprintf (stderr, "WARNING: NAN in vector1\n");
     833            }
     834        }
    779835        {
    780836            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
    … …  
    798854
    799855        psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    800         psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
     856        psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    801857        psFree(sumMatrix);
    802858        if (!luMatrix) {
    … …  
    807863            return NULL;
    808864        }
    809         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
     865        kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
     866
     867#ifdef TESTING
     868        // XXX double-check for NAN in data:
     869        for (int ix = 0; ix < kernels->solution1->n; ix++) {
     870            if (!isfinite(kernels->solution1->data.F64[ix])) {
     871                fprintf (stderr, "WARNING: NAN in vector1\n");
     872            }
     873        }
     874#endif
    810875
    811876        psFree(sumVector);
    … …  
    10621127                source = stamp->image1;
    10631128                convolutions = stamp->convolutions1;
     1129
     1130                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
     1131                // so that it is on the scale of image1.
     1132                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
     1133                                                          false); // Kernel image
     1134                if (!image) {
     1135                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     1136                    return false;
     1137                }
     1138                double sumKernel = 0;   // Sum of kernel, for normalising residual
     1139                int size = kernels->size; // Half-size of kernel
     1140                int fullSize = 2 * size + 1; // Full size of kernel
     1141                for (int y = 0; y < fullSize; y++) {
     1142                    for (int x = 0; x < fullSize; x++) {
     1143                        sumKernel += image->data.F32[y][x];
     1144                    }
     1145                }
     1146                psFree(image);
     1147                devNorm = 1.0 / sumKernel / numPixels;
    10641148                break;
    10651149              case PM_SUBTRACTION_MODE_2:

  • branches/pap/psModules/src/imcombine/pmSubtractionKernels.c

    r20421 r25027  
    145145
    146146                // Calculate moments
    147                 double sum = 0.0;       // Sum of kernel component, for normalisation
    148147                double moment = 0.0;    // Moment, for penalty
    149148                for (int v = -size, y = 0; v <= size; v++, y++) {
    150149                    for (int u = -size, x = 0; u <= size; u++, x++) {
    151150                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
    152                         sum += value;
    153151                        moment += value * (PS_SQR(u) + PS_SQR(v));
    154152                    }
    … …  
    163161                        }
    164162                    }
    165                     sum = 1.0 / sum;
     163                    sum = 1.0 / sqrt(sum);
    166164                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    167165                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    168                     moment *= sum;
     166                    moment *= PS_SQR(sum);
    169167                }
    170168

  • branches/pap/psModules/src/imcombine/pmSubtractionMask.c

    r23851 r25027  
    66#include <pslib.h>
    77
     8#include "pmErrorCodes.h"
    89#include "pmSubtraction.h"
    910#include "pmSubtractionKernels.h"
    … …  
    7879        }
    7980        if (numBad > badFrac * numCols * numRows) {
    80             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     81            psError(PM_ERR_SMALL_AREA, true,
    8182                    "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
    8283                    numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac);

  • branches/pap/psModules/src/imcombine/pmSubtractionMatch.c

    r23944 r25027  
    246246                                         badFrac, mode); // Subtraction mask
    247247    if (!subMask) {
    248         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
     248        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    249249        psFree(kernels);
    250250        psFree(regions);
    … …  
    337337        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false);
    338338    }
    339     PS_ASSERT_INT_POSITIVE(iter, false);
     339    PS_ASSERT_INT_NONNEGATIVE(iter, false);
    340340    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    341341
    … …  
    379379                                badFrac, subMode);
    380380    if (!subMask) {
    381         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
     381        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    382382        goto MATCH_ERROR;
    383383    }
    … …  
    432432            }
    433433
    434             if (sources) {
     434            if (stampsName && strlen(stampsName) > 0) {
     435                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
     436                                                        footprint, stampSpacing, subMode);
     437            } else if (sources) {
    435438                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    436439                                                           footprint, stampSpacing, subMode);
    437             } else if (stampsName && strlen(stampsName) > 0) {
    438                 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    439                                                         footprint, stampSpacing, subMode);
    440440            }
    441441
    … …  
    830830    psFree(mask);
    831831
     832    // XXX raise an error here or not?
     833    if (isnan(stats->robustMedian)) {
     834        psFree(stats);
     835        return PM_SUBTRACTION_MODE_ERR;
     836    }
     837
    832838    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median width ratio: %lf", stats->robustMedian);
    833839    pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
    834840    psFree(stats);
    835841
    836     // XXX EAM : I think Paul left some test code in here.  I've commented these lines out
    837     // return PM_SUBTRACTION_MODE_2;
    838     // exit(1);
    839 
    840842    return mode;
    841843}

  • branches/pap/psModules/src/imcombine/pmSubtractionStamps.c

    r23937 r25027  
    7878}
    7979
    80 // Is this position unmasked?
    81 static bool checkStampMask(int x, int y, // Coordinates of stamp
    82                            int numCols, int numRows, // Size of image
    83                            const psImage *mask, // Mask
    84                            pmSubtractionMode mode, // Mode for subtraction
    85                            int footprint, // Size of stamp
    86                            int border // Size of border
    87                            )
    88 {
    89     // Check the footprint bounds
    90     if (x < border || x >= numCols - border || y < border || y >= numRows - border) {
    91         return false;
    92     }
    93 
    94     if (!mask) {
    95         return true;
    96     }
    97 
    98     // Check the central pixel
    99     if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
    100         mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
    101         return false;
    102     }
    103 
    104     return true;
     80
     81// Search a region for a suitable stamp
     82bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return
     83                 float *fluxStamp, // Flux of stamp, to return
     84                 const psImage *image1, const psImage *image2, // Images to search
     85                 float thresh1, float thresh2, // Thresholds for images
     86                 const psImage *subMask, // Subtraction mask
     87                 int xMin, int xMax, int yMin, int yMax, // Bounds of search
     88                 int numCols, int numRows, // Size of images
     89                 int border             // Border around image
     90    )
     91{
     92    bool found = false;                 // Found a suitable stamp?
     93    *fluxStamp = -INFINITY;             // Flux of best stamp
     94
     95    // Ensure we're not going to go outside the bounds of the image
     96    xMin = PS_MAX(border, xMin);
     97    xMax = PS_MIN(numCols - border - 1, xMax);
     98    yMin = PS_MAX(border, yMin);
     99    yMax = PS_MIN(numRows - border - 1, yMax);
     100
     101    for (int y = yMin; y <= yMax; y++) {
     102        for (int x = xMin; x <= xMax; x++) {
     103            if ((image1 && image1->data.F32[y][x] < thresh1) ||
     104                (image2 && image2->data.F32[y][x] < thresh2)) {
     105                continue;
     106            }
     107
     108            if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] &
     109                (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
     110                return false;
     111            }
     112
     113            // We take the MIN to attempt to avoid transients in both images
     114            float flux = (image1 && image2) ? PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]) :
     115                ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel
     116            if (flux > *fluxStamp) {
     117                *xStamp = x;
     118                *yStamp = y;
     119                *fluxStamp = flux;
     120                found = true;
     121            }
     122        }
     123    }
     124
     125    return found;
    105126}
    106127
    … …  
    314335            numSearch++;
    315336
    316             float xStamp = 0, yStamp = 0;  // Coordinates of stamp
    317             float fluxStamp = NAN;          // Flux of stamp
    318             bool goodStamp = false;         // Found a good stamp?
     337            int xStamp = 0, yStamp = 0; // Coordinates of stamp
     338            float fluxStamp = -INFINITY; // Flux of stamp
     339            bool goodStamp = false;     // Found a good stamp?
    319340
    320341            // A couple different ways of finding stamps:
    … …  
    324345                psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
    325346
    326                 // Take stamp off the top of the (sorted) list
    327                 if (xList->n > 0) {
    328                     int index = xList->n - 1; // Index of new stamp
    329                     xStamp = xList->data.F32[index];
    330                     yStamp = yList->data.F32[index];
    331                     fluxStamp = fluxList->data.F32[index];
     347                // Take stamps off the top of the (sorted) list
     348                for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
     349                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
    332350
    333351                    // Chop off the top of the list
    334                     xList->n = index;
    335                     yList->n = index;
    336                     fluxList->n = index;
    337 
    338                     goodStamp = true;
    339                     int x = xStamp + 0.5, y = yStamp + 0.5; // Pixel location
    340                     if (image1 && image1->data.F32[y][x] < thresh1) {
    341                         goodStamp = false;
    342                     } else if (image2 && image2->data.F32[y][x] < thresh2) {
    343                         goodStamp = false;
    344                     }
    345                 } else {
    346                     psTrace("psModules.imcombine", 9, "No sources in subregion %d", i);
     352                    xList->n = j;
     353                    yList->n = j;
     354                    fluxList->n = j;
     355
     356                    // Fish around a bit to see if we can find a pixel that isn't masked
     357                    psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
     358                            i, xCentre, yCentre);
     359
     360                    // Search bounds
     361                    int search = footprint - size; // Search radius
     362                    int xMin = PS_MAX(border, xCentre - search);
     363                    int xMax = PS_MIN(numCols - border -1, xCentre + search);
     364                    int yMin = PS_MAX(border, yCentre - search);
     365                    int yMax = PS_MIN(numRows - border - 1, yCentre + search);
     366
     367                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
     368                                            subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
    347369                }
    348370            } else {
    349371                // Use a simple method of automatically finding stamps --- take the highest pixel in the
    350372                // subregion
    351                 fluxStamp = 0.0;
    352373                psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
    353374
    354                 const psImage *image = NULL; // Image for flux of stamp
    355                 switch (mode) {
    356                   case PM_SUBTRACTION_MODE_1:
    357                     image = image2 ? image2 : image1;
    358                     break;
    359                   case PM_SUBTRACTION_MODE_2:
    360                     image = image1 ? image1 : image2;
    361                     break;
    362                   case PM_SUBTRACTION_MODE_UNSURE:
    363                   case PM_SUBTRACTION_MODE_DUAL:
    364                     // If both are available, we'll take the MIN value below in the hope of rejecting
    365                     // transients in both
    366                     image = (image1 && image2) ? NULL : (image1 ? image1 : image2);
    367                   default:
    368                     psAbort("Unrecognised mode: %x", mode);
    369                 }
    370 
    371                 for (int y = subRegion->y0; y <= subRegion->y1; y++) {
    372                     for (int x = subRegion->x0; x <= subRegion->x1; x++) {
    373                         if ((image1 && image1->data.F32[y][x] < thresh1) ||
    374                             (image2 && image2->data.F32[y][x] < thresh2)) {
    375                             continue;
    376                         }
    377                         float value = image ? image->data.F32[y][x] :
    378                             PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]); // Value of image
    379                         if (value > fluxStamp &&
    380                             checkStampMask(x, y, numCols, numRows, subMask, mode, footprint, border)) {
    381                             fluxStamp = image->data.F32[y][x];
    382                             xStamp = x;
    383                             yStamp = y;
    384                             goodStamp = true;
    385                         }
    386                     }
    387                 }
     375                goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
     376                                        subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1,
     377                                        numCols, numRows, border);
    388378            }
    389379
    … …  
    453443                                                                 region, footprint, spacing); // Stamp list
    454444    int numStamps = stamps->num;        // Number of stamps
    455     int numCols = image->numCols, numRows = image->numRows; // Size of image
    456     int border = footprint + size;      // Border size
    457445
    458446    psString ds9name = NULL;            // Filename for ds9 region file
    … …  
    480468            psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region",
    481469                    xPix, yPix);
    482             pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta");
    483             continue;
    484         }
    485         if (!checkStampMask(xPix, yPix, numCols, numRows, subMask, mode, footprint, border)) {
    486             // Not a good stamp
    487             psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",
    488                     xPix, yPix);
    489470            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red");
    490471            continue;

  • branches/pap/psModules/src/objects/Makefile.am

    r23184 r25027  
    1010        pmFootprintArraysMerge.c \
    1111        pmFootprintAssignPeaks.c \
     12        pmFootprintCullPeaks.c \
    1213        pmFootprintFind.c \
    1314        pmFootprintFindAtPoint.c \
    … …  
    3839        pmSourceIO_PS1_CAL_0.c \
    3940        pmSourceIO_CMF_PS1_V1.c \
     41        pmSourceIO_CMF_PS1_V2.c \
     42        pmSourceIO_MatchedRefs.c \
    4043        pmSourcePlots.c \
    4144        pmSourcePlotPSFModel.c \

  • branches/pap/psModules/src/objects/models/pmModel_PS1_V1.c

    r20001 r25027  
    159159            break;
    160160        case PM_PAR_7:
    161             params_min =   0.1;
     161            params_min =  -1.0;
    162162            break;
    163163        default:

  • branches/pap/psModules/src/objects/models/pmModel_SGAUSS.c

    r15834 r25027  
    325325
    326326    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    327     psVectorStats (stats, contour, NULL, NULL, 0);
     327    if (!psVectorStats (stats, contour, NULL, NULL, 0)) {
     328        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     329        return false;
     330    }
    328331    value = stats->sampleMedian;
    329332

  • branches/pap/psModules/src/objects/pmFootprint.h

    r20945 r25027  
    5151void pmSetFootprintArrayIDsForImage(psImage *idImage,
    5252                                    const psArray *footprints, // the footprints to insert
    53                                     const bool relativeIDs); // show IDs starting at 0, not pmFootprint->id
     53                                    const bool relativeIDs // show IDs starting at 0, not pmFootprint->id
     54    );
    5455
    5556psErrorCode pmFootprintsAssignPeaks(psArray *footprints, const psArray *peaks);
    5657
    57 psErrorCode pmFootprintArrayCullPeaks(const psImage *img, const psImage *weight, psArray *footprints,
    58                                  const float nsigma, const float threshold_min);
    59 psErrorCode pmFootprintCullPeaks(const psImage *img, const psImage *weight, pmFootprint *fp,
    60                                  const float nsigma, const float threshold_min);
     58psErrorCode pmFootprintCullPeaks(const psImage *img,       // the image wherein lives the footprint
     59                                 const psImage *weight,    // corresponding variance image
     60                                 pmFootprint *fp,          // Footprint containing mortal peaks
     61                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
     62                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
     63                                 const float min_threshold // minimum permitted coll height
     64    );
    6165
    6266psArray *pmFootprintArrayToPeaks(const psArray *footprints);

  • branches/pap/psModules/src/objects/pmFootprintArraysMerge.c

    r20937 r25027  
    2828 */
    2929psArray *pmFootprintArraysMerge(const psArray *footprints1, // one set of footprints
    30                                 const psArray *footprints2, // the other set
    31                                 const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2
    32     assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0]));
    33     assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0]));
     30                                const psArray *footprints2, // the other set
     31                                const int includePeaks // which peaks to set? 0x1 => footprints1, 0x2 => 2
     32    )
     33{
     34    if (!footprints1 && !footprints2) {
     35        // No footprints in merged array
     36        return psArrayAllocEmpty(0);
     37    }
    3438
    35     if (footprints1->n == 0 || footprints2->n == 0) {           // nothing to do but put copies on merged
    36         const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1;
     39    assert(!footprints1 || footprints1->n == 0 || pmFootprintTest(footprints1->data[0]));
     40    assert(!footprints2 || footprints2->n == 0 || pmFootprintTest(footprints2->data[0]));
    3741
    38         psArray *merged = psArrayAllocEmpty(old->n);
    39         for (int i = 0; i < old->n; i++) {
    40             psArrayAdd(merged, 1, old->data[i]);
    41         }
    42        
    43         return merged;
     42    if (!footprints1 || footprints1->n == 0 || !footprints2 || footprints2->n == 0) {
     43        // nothing to do but put copies on merged
     44        const psArray *old = (!footprints1 || footprints1->n == 0) ? footprints2 : footprints1;
     45
     46        psArray *merged = psArrayAllocEmpty(old->n);
     47        for (int i = 0; i < old->n; i++) {
     48            psArrayAdd(merged, 1, old->data[i]);
     49        }
     50
     51        return merged;
    4452    }
    4553    /*
    … …  
    4856     */
    4957    {
    50         pmFootprint *fp1 = footprints1->data[0];
    51         pmFootprint *fp2 = footprints2->data[0];
    52         if (fp1->region.x0 != fp2->region.x0 ||
    53             fp1->region.x1 != fp2->region.x1 ||
    54             fp1->region.y0 != fp2->region.y0 ||
    55             fp1->region.y1 != fp2->region.y1) {
    56             psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    57                     "The two pmFootprint arrays correspnond to different-sized regions");
    58             return NULL;
    59         }
     58        pmFootprint *fp1 = footprints1->data[0];
     59        pmFootprint *fp2 = footprints2->data[0];
     60        if (fp1->region.x0 != fp2->region.x0 ||
     61            fp1->region.x1 != fp2->region.x1 ||
     62            fp1->region.y0 != fp2->region.y0 ||
     63            fp1->region.y1 != fp2->region.y1) {
     64            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     65                    "The two pmFootprint arrays correspnond to different-sized regions");
     66            return NULL;
     67        }
    6068    }
    6169    /*
    … …  
    7280     * Now assign the peaks appropriately.  We could do this more efficiently
    7381     * using idImage (which we just freed), but this is easy and probably fast enough
    74      */ 
     82     */
    7583    if (includePeaks & 0x1) {
    76         psArray *peaks = pmFootprintArrayToPeaks(footprints1);
    77         pmFootprintsAssignPeaks(merged, peaks);
    78         psFree(peaks);
     84        psArray *peaks = pmFootprintArrayToPeaks(footprints1);
     85        pmFootprintsAssignPeaks(merged, peaks);
     86        psFree(peaks);
    7987    }
    8088
    8189    if (includePeaks & 0x2) {
    82         psArray *peaks = pmFootprintArrayToPeaks(footprints2);
    83         pmFootprintsAssignPeaks(merged, peaks);
    84         psFree(peaks);
     90        psArray *peaks = pmFootprintArrayToPeaks(footprints2);
     91        pmFootprintsAssignPeaks(merged, peaks);
     92        psFree(peaks);
    8593    }
    86    
     94
    8795    return merged;
    8896}

  • branches/pap/psModules/src/objects/pmGrowthCurveGenerate.c

    r21183 r25027  
    8989        psVectorAppend (values, growth->fitMag);
    9090    }
    91     psVectorStats (stats, values, NULL, NULL, 0);
     91    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
     92        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     93        return false;
     94    }
    9295    psf->growth->fitMag = stats->sampleMedian;
    9396
    … …  
    104107            psVectorAppend (values, growth->apMag->data.F32[i]);
    105108        }
    106         psVectorStats (stats, values, NULL, NULL, 0);
     109        if (!psVectorStats (stats, values, NULL, NULL, 0)) {
     110            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     111            return false;
     112        }
    107113        psf->growth->apMag->data.F32[i] = stats->sampleMedian;
    108114    }

  • branches/pap/psModules/src/objects/pmMoments.c

    r23487 r25027  
    5151    tmp->Sky = 0.0;
    5252    tmp->nPixels = 0;
     53    tmp->SN = 0;
    5354
    5455    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);

  • branches/pap/psModules/src/objects/pmPSF.c

    r23487 r25027  
    7575    options->poissonErrorsParams  = true;
    7676
     77    options->chiFluxTrend = true;
     78
    7779    return options;
    7880}

  • branches/pap/psModules/src/objects/pmPSF.h

    r23487 r25027  
    8282    bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
    8383    float         radius;
     84    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
    8485} pmPSFOptions;
    8586

  • branches/pap/psModules/src/objects/pmPSFtry.c

    r21183 r25027  
    4343
    4444    for (int j = 0; j < trend->map->map->numRows; j++) {
    45         for (int i = 0; i < trend->map->map->numCols; i++) {
    46             fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
    47         }
    48         fprintf (stderr, "\t\t\t");
    49         for (int i = 0; i < trend->map->map->numCols; i++) {
    50             fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
    51         }
    52         fprintf (stderr, "\n");
     45        for (int i = 0; i < trend->map->map->numCols; i++) {
     46            fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
     47        }
     48        fprintf (stderr, "\t\t\t");
     49        for (int i = 0; i < trend->map->map->numCols; i++) {
     50            fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
     51        }
     52        fprintf (stderr, "\n");
    5353    }
    5454    return true;
    … …  
    6363    float Wt = 0.0;
    6464    for (int j = 0; j < map->map->numRows; j++) {
    65         for (int i = 0; i < map->map->numCols; i++) {
    66             if (!isfinite(map->map->data.F32[j][i])) continue;
    67             Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
    68             Wt += map->error->data.F32[j][i];
    69         }
     65        for (int i = 0; i < map->map->numCols; i++) {
     66            if (!isfinite(map->map->data.F32[j][i])) continue;
     67            Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
     68            Wt += map->error->data.F32[j][i];
     69        }
    7070    }
    7171
    … …  
    7575    // XXX for now, we are just replacing bad pixels with the Mean
    7676    for (int j = 0; j < map->map->numRows; j++) {
    77         for (int i = 0; i < map->map->numCols; i++) {
    78             if (isfinite(map->map->data.F32[j][i]) &&
    79                 (map->error->data.F32[j][i] > 0.0)) continue;
    80             map->map->data.F32[j][i] = Mean;
    81         }
     77        for (int i = 0; i < map->map->numCols; i++) {
     78            if (isfinite(map->map->data.F32[j][i]) &&
     79                (map->error->data.F32[j][i] > 0.0)) continue;
     80            map->map->data.F32[j][i] = Mean;
     81        }
    8282    }
    8383    return true;
    … …  
    174174
    175175        pmSource *source = psfTry->sources->data[i];
    176         if (!source->moments) {
     176        if (!source->moments) {
    177177            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    178             continue;
    179         }
    180         if (!source->moments->nPixels) {
     178            continue;
     179        }
     180        if (!source->moments->nPixels) {
    181181            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    182             continue;
    183         }
     182            continue;
     183        }
    184184
    185185        source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
    … …  
    211211
    212212    if (Next == 0) {
    213         psError(PS_ERR_UNKNOWN, true, "No sources with good extended fits from which to determine PSF.");
     213        psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF.");
    214214        psFree(psfTry);
    215215        return NULL;
    … …  
    282282
    283283    if (Npsf == 0) {
    284         psError(PS_ERR_UNKNOWN, true, "No sources with good PSF fits after model is built.");
     284        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
    285285        psFree(psfTry);
    286286        return NULL;
    … …  
    311311
    312312    // linear clipped fit of chisq trend vs flux
    313     bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, mask, 0xff, chisq, NULL, flux);
    314     psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    315     psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    316 
    317     psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
    318               psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
    319 
    320     psFree(flux);
    321     psFree(mask);
    322     psFree(chisq);
    323 
    324     if (!result) {
    325         psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
    326         psFree(psfTry);
    327         return NULL;
     313    if (options->chiFluxTrend) {
     314        bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
     315                                                  mask, 0xff, chisq, NULL, flux);
     316        psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
     317        psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
     318
     319        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
     320                  psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
     321
     322        psFree(flux);
     323        psFree(mask);
     324        psFree(chisq);
     325
     326        if (!result) {
     327            psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
     328            psFree(psfTry);
     329            return NULL;
     330        }
    328331    }
    329332
    … …  
    600603
    601604    for (int i = 0; i < sources->n; i++) {
    602         // skip any masked sources (failed to fit one of the model steps or get a magnitude)
    603         if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    604        
     605        // skip any masked sources (failed to fit one of the model steps or get a magnitude)
     606        if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     607
    605608        pmSource *source = sources->data[i];
    606609        assert (source->modelEXT); // all unmasked sources should have modelEXT
    … …  
    628631        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    629632
    630             // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    631             for (int i = 0; i < mask->n; i++) {
    632                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    633             }
     633            // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
     634            for (int i = 0; i < mask->n; i++) {
     635                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     636            }
    634637            if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    635638                break;
    … …  
    643646        }
    644647        if (entryMin == -1) {
    645             psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params");
     648            psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params");
    646649            return false;
    647650        }
    648651
    649         // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    650         for (int i = 0; i < mask->n; i++) {
    651             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    652         }
     652        // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
     653        for (int i = 0; i < mask->n; i++) {
     654            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     655        }
    653656        if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    654657            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    … …  
    659662        psf->trendNy = trend->map->map->numRows;
    660663
    661         // copy mask back to srcMask
    662         for (int i = 0; i < mask->n; i++) {
    663             srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    664         }
     664        // copy mask back to srcMask
     665        for (int i = 0; i < mask->n; i++) {
     666            srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     667        }
    665668
    666669        psFree (mask);
    … …  
    686689        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    687690
    688             psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    689             psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    690 
    691             pmTrend2D *trend = NULL;
    692             float mean, stdev;
     691            psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     692            psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     693
     694            pmTrend2D *trend = NULL;
     695            float mean, stdev;
    693696
    694697            // XXX we are using the same stats structure on each pass: do we need to re-init it?
    695698            bool status = true;
    696699
    697             trend = psf->params->data[PM_PAR_E0];
     700            trend = psf->params->data[PM_PAR_E0];
    698701            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
    699             mean = psStatsGetValue (trend->stats, meanOption);
    700             stdev = psStatsGetValue (trend->stats, stdevOption);
     702            mean = psStatsGetValue (trend->stats, meanOption);
     703            stdev = psStatsGetValue (trend->stats, stdevOption);
    701704            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    702             pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    703 
    704             trend = psf->params->data[PM_PAR_E1];
     705            pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
     706
     707            trend = psf->params->data[PM_PAR_E1];
    705708            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
    706             mean = psStatsGetValue (trend->stats, meanOption);
    707             stdev = psStatsGetValue (trend->stats, stdevOption);
     709            mean = psStatsGetValue (trend->stats, meanOption);
     710            stdev = psStatsGetValue (trend->stats, stdevOption);
    708711            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    709             pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    710 
    711             trend = psf->params->data[PM_PAR_E2];
     712            pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
     713
     714            trend = psf->params->data[PM_PAR_E2];
    712715            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
    713             mean = psStatsGetValue (trend->stats, meanOption);
    714             stdev = psStatsGetValue (trend->stats, stdevOption);
     716            mean = psStatsGetValue (trend->stats, meanOption);
     717            stdev = psStatsGetValue (trend->stats, stdevOption);
    715718            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    716             pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
     719            pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    717720
    718721            if (!status) {
    … …  
    755758    if (psf->fieldNx > psf->fieldNy) {
    756759        Nx = scale;
    757         float AR = psf->fieldNy / (float) psf->fieldNx;
    758         Ny = (int) (Nx * AR + 0.5);
     760        float AR = psf->fieldNy / (float) psf->fieldNx;
     761        Ny = (int) (Nx * AR + 0.5);
    759762        Ny = PS_MAX (1, Ny);
    760763    } else {
    761764        Ny = scale;
    762         float AR = psf->fieldNx / (float) psf->fieldNy;
    763         Nx = (int) (Ny * AR + 0.5);
     765        float AR = psf->fieldNx / (float) psf->fieldNy;
     766        Nx = (int) (Ny * AR + 0.5);
    764767        Nx = PS_MAX (1, Nx);
    765768    }
    … …  
    809812
    810813    if (scale == 1) {
    811         x_fit  = psMemIncrRefCounter (x);
    812         y_fit  = psMemIncrRefCounter (y);
    813         x_tst  = psMemIncrRefCounter (x);
    814         y_tst  = psMemIncrRefCounter (y);
    815         e0obs_fit = psMemIncrRefCounter (e0obs);
    816         e1obs_fit = psMemIncrRefCounter (e1obs);
    817         e2obs_fit = psMemIncrRefCounter (e2obs);
    818         e0obs_tst = psMemIncrRefCounter (e0obs);
    819         e1obs_tst = psMemIncrRefCounter (e1obs);
    820         e2obs_tst = psMemIncrRefCounter (e2obs);
     814        x_fit  = psMemIncrRefCounter (x);
     815        y_fit  = psMemIncrRefCounter (y);
     816        x_tst  = psMemIncrRefCounter (x);
     817        y_tst  = psMemIncrRefCounter (y);
     818        e0obs_fit = psMemIncrRefCounter (e0obs);
     819        e1obs_fit = psMemIncrRefCounter (e1obs);
     820        e2obs_fit = psMemIncrRefCounter (e2obs);
     821        e0obs_tst = psMemIncrRefCounter (e0obs);
     822        e1obs_tst = psMemIncrRefCounter (e1obs);
     823        e2obs_tst = psMemIncrRefCounter (e2obs);
    821824    } else {
    822         x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    823         y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    824         x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    825         y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    826         e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    827         e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    828         e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    829         e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    830         e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    831         e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    832         for (int i = 0; i < e0obs_fit->n; i++) {
    833             // e0obs->n ==  8 or 9:
    834             // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
    835             // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
    836             x_fit->data.F32[i] = x->data.F32[2*i+0]; 
    837             x_tst->data.F32[i] = x->data.F32[2*i+1]; 
    838             y_fit->data.F32[i] = y->data.F32[2*i+0]; 
    839             y_tst->data.F32[i] = y->data.F32[2*i+1]; 
    840 
    841             e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; 
    842             e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; 
    843             e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
    844             e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
    845             e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
    846             e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
    847         }
     825        x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     826        y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     827        x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     828        y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     829        e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     830        e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     831        e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     832        e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     833        e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     834        e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     835        for (int i = 0; i < e0obs_fit->n; i++) {
     836            // e0obs->n ==  8 or 9:
     837            // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
     838            // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
     839            x_fit->data.F32[i] = x->data.F32[2*i+0];
     840            x_tst->data.F32[i] = x->data.F32[2*i+1];
     841            y_fit->data.F32[i] = y->data.F32[2*i+0];
     842            y_tst->data.F32[i] = y->data.F32[2*i+1];
     843
     844            e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];
     845            e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];
     846            e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
     847            e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
     848            e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
     849            e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
     850        }
    848851    }
    849852
    … …  
    852855    // copy mask values to fitMask as a starting point
    853856    for (int i = 0; i < fitMask->n; i++) {
    854         fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     857        fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    855858    }
    856859
    … …  
    859862    for (int i = 0; i < nIter; i++) {
    860863        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    861         psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    862         psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    863 
    864         pmTrend2D *trend = NULL;
    865         float mean, stdev;
    866 
    867         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    868         bool status = true;
    869 
    870         trend = psf->params->data[PM_PAR_E0];
    871         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    872         mean = psStatsGetValue (trend->stats, meanOption);
    873         stdev = psStatsGetValue (trend->stats, stdevOption);
    874         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    875         // printTrendMap (trend);
    876         psImageMapCleanup (trend->map);
    877         // printTrendMap (trend);
    878         pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    879 
    880         trend = psf->params->data[PM_PAR_E1];
    881         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    882         mean = psStatsGetValue (trend->stats, meanOption);
    883         stdev = psStatsGetValue (trend->stats, stdevOption);
    884         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    885         // printTrendMap (trend);
    886         psImageMapCleanup (trend->map);
    887         // printTrendMap (trend);
    888         pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    889 
    890         trend = psf->params->data[PM_PAR_E2];
    891         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    892         mean = psStatsGetValue (trend->stats, meanOption);
    893         stdev = psStatsGetValue (trend->stats, stdevOption);
    894         psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    895         // printTrendMap (trend);
    896         psImageMapCleanup (trend->map);
    897         // printTrendMap (trend);
    898         pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
     864        psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     865        psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     866
     867        pmTrend2D *trend = NULL;
     868        float mean, stdev;
     869
     870        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     871        bool status = true;
     872
     873        trend = psf->params->data[PM_PAR_E0];
     874        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
     875        mean = psStatsGetValue (trend->stats, meanOption);
     876        stdev = psStatsGetValue (trend->stats, stdevOption);
     877        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
     878        // printTrendMap (trend);
     879        psImageMapCleanup (trend->map);
     880        // printTrendMap (trend);
     881        pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
     882
     883        trend = psf->params->data[PM_PAR_E1];
     884        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
     885        mean = psStatsGetValue (trend->stats, meanOption);
     886        stdev = psStatsGetValue (trend->stats, stdevOption);
     887        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
     888        // printTrendMap (trend);
     889        psImageMapCleanup (trend->map);
     890        // printTrendMap (trend);
     891        pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
     892
     893        trend = psf->params->data[PM_PAR_E2];
     894        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
     895        mean = psStatsGetValue (trend->stats, meanOption);
     896        stdev = psStatsGetValue (trend->stats, stdevOption);
     897        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
     898        // printTrendMap (trend);
     899        psImageMapCleanup (trend->map);
     900        // printTrendMap (trend);
     901        pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    899902    }
    900903    psf->psfTrendStats->clipIter = nIter; // restore default setting
    … …  
    941944    // XXX copy fitMask values back to mask
    942945    for (int i = 0; i < fitMask->n; i++) {
    943         mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     946        mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    944947    }
    945948    psFree (fitMask);
    … …  
    958961    float dEsquare = 0.0;
    959962    psStatsInit (stats);
    960     psVectorStats (stats, e0res, NULL, mask, maskValue);
     963    if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
     964        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     965        return false;
     966    }
    961967    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    962968
    963969    psStatsInit (stats);
    964     psVectorStats (stats, e1res, NULL, mask, maskValue);
     970    if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
     971        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     972        return false;
     973    }
    965974    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    966975
    967976    psStatsInit (stats);
    968     psVectorStats (stats, e2res, NULL, mask, maskValue);
     977    if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
     978        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     979        return false;
     980    }
    969981    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    970982
    … …  
    9981010    for (int i = 0; i < nBin; i++) {
    9991011        int j;
    1000         int nValid = 0;
     1012        int nValid = 0;
    10011013        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    10021014            int N = index->data.U32[n];
    … …  
    10061018
    10071019            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    1008             if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
    1009         }
    1010         if (nValid < 3) continue;
     1020            if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
     1021        }
     1022        if (nValid < 3) continue;
    10111023
    10121024        dE0subset->n = j;
    … …  
    10181030        float dEsquare = 0.0;
    10191031        psStatsInit (statsS);
    1020         psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
     1032        if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
     1033        }
    10211034        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10221035
    10231036        psStatsInit (statsS);
    1024         psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
     1037        if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
     1038            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1039            return false;
     1040        }
    10251041        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10261042
    10271043        psStatsInit (statsS);
    1028         psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
     1044        if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
     1045            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1046            return false;
     1047        }
    10291048        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10301049

  • branches/pap/psModules/src/objects/pmPeaks.c

    r23187 r25027  
    147147    tmp->y = y;
    148148    tmp->value = value;
    149     tmp->flux = 0;
     149    tmp->flux = value;
    150150    tmp->SN = 0;
    151151    tmp->xf = x;

  • branches/pap/psModules/src/objects/pmSource.c

    r23187 r25027  
    242242    extend |= (mySource->maskView == NULL);
    243243
    244     // extend = true;
    245244    if (extend) {
    246245        // re-create the subimage
    … …  
    250249
    251250        mySource->pixels   = psImageSubset(readout->image,  newRegion);
    252         mySource->variance   = psImageSubset(readout->variance, newRegion);
     251        mySource->variance = psImageSubset(readout->variance, newRegion);
    253252        mySource->maskView = psImageSubset(readout->mask,   newRegion);
    254253        mySource->region   = newRegion;
    … …  
    291290
    292291    bool status = true;                 // Status of MD lookup
    293     float PSF_CLUMP_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_SN_LIM");
     292    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    294293    if (!status) {
    295         PSF_CLUMP_SN_LIM = 0;
     294        PSF_SN_LIM = 0;
    296295    }
    297296    float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
    … …  
    342341            }
    343342
    344             if (src->moments->SN < PSF_CLUMP_SN_LIM) {
     343            if (src->moments->SN < PSF_SN_LIM) {
    345344                psTrace("psModules.objects", 10, "Rejecting source from clump because of low S/N (%f)\n",
    346345                        src->moments->SN);
    … …  
    450449            if (tmpSrc->moments == NULL)
    451450                continue;
    452             if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM)
     451            if (tmpSrc->moments->SN < PSF_SN_LIM)
    453452                continue;
    454453
    … …  
    479478        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    480479
    481         psVectorStats (stats, tmpSx, NULL, NULL, 0);
     480        if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) {
     481            psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
     482            return (emptyClump);
     483        }
    482484        psfClump.X  = stats->clippedMean;
    483485        psfClump.dX = stats->clippedStdev;
    484486
    485         psVectorStats (stats, tmpSy, NULL, NULL, 0);
     487        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
     488            psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
     489            return (emptyClump);
     490        }
    486491        psfClump.Y  = stats->clippedMean;
    487492        psfClump.dY = stats->clippedStdev;
    … …  
    528533    bool status;
    529534    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
    530     if (!status)
    531         PSF_SN_LIM = 20.0;
     535    if (!status) PSF_SN_LIM = 20.0;
    532536    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    533     if (!status)
    534         PSF_CLUMP_NSIGMA = 1.5;
    535     float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     537    if (!status) PSF_CLUMP_NSIGMA = 1.5;
     538
     539    // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    536540
    537541    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    … …  
    543547        pmSource *source = (pmSource *) sources->data[i];
    544548
    545         // psf clumps are found for image subregions:
    546         // skip sources not in this region
     549        // psf clumps are found for image subregions:
     550        // skip sources not in this region
    547551        if (source->peak->x <  region->x0) continue;
    548552        if (source->peak->x >= region->x1) continue;
    549553        if (source->peak->y <  region->y0) continue;
    550         if (source->peak->y >= region->y1) continue;
    551 
    552         // should be set by pmSourceAlloc
     554        if (source->peak->y >= region->y1) continue;
     555
     556        // should be set by pmSourceAlloc
    553557        psAssert (source->type == PM_SOURCE_TYPE_UNKNOWN, "source type was not init-ed?");
    554558
    … …  
    556560        if (!source->moments) {
    557561            source->type = PM_SOURCE_TYPE_STAR;
    558             psAssert (source->mode & noMoments, "why is this source missing moments?");
     562            psAssert (source->mode & noMoments, "why is this source missing moments?");
    559563            Nstar++;
    560564            continue;
    … …  
    576580            source->type = PM_SOURCE_TYPE_STAR;
    577581            source->mode |= PM_SOURCE_MODE_SATSTAR;
    578             // recalculate moments here with larger box?
    579             pmSourceMoments (source, INNER_RADIUS);
     582            // why do we recalculate moments here?
     583            // we already attempt to do this in psphotSourceStats
     584            // pmSourceMoments (source, INNER_RADIUS);
    580585            Nsatstar ++;
    581586            continue;
    … …  
    590595        }
    591596
    592         // likely defect (too small to be stellar) (push out to 3 sigma)
    593         // low S/N objects which are small are probably stellar
    594         // only set candidate defects if
    595         // XXX these limits are quite arbitrary
    596         if ((sigX < 0.05) || (sigY < 0.05)) {
    597             source->type = PM_SOURCE_TYPE_DEFECT;
    598             source->mode |= PM_SOURCE_MODE_DEFECT;
    599             Ncr ++;
    600             continue;
    601         }
    602 
    603         // likely unsaturated extended source (too large to be stellar)
    604         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    605             source->type = PM_SOURCE_TYPE_EXTENDED;
    606             Next ++;
    607             continue;
    608         }
    609 
    610         // the rest are probable stellar objects
    611         starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
    612         starsn_moments->n ++;
    613         starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
    614         starsn_peaks->n ++;
    615         Nstar ++;
    616 
    617         // PSF star (within 1.5 sigma of clump center, S/N > limit)
    618         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    619         if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
    620             source->type = PM_SOURCE_TYPE_STAR;
    621             source->mode |= PM_SOURCE_MODE_PSFSTAR;
    622             Npsf ++;
    623             continue;
     597        // The following determinations require the use of moments
     598        if (!(source->mode & noMoments)) {
     599            // likely defect (too small to be stellar) (push out to 3 sigma)
     600            // low S/N objects which are small are probably stellar
     601            // XXX these limits are quite arbitrary
     602            if (sigX < 0.05 || sigY < 0.05) {
     603                source->type = PM_SOURCE_TYPE_DEFECT;
     604                source->mode |= PM_SOURCE_MODE_DEFECT;
     605                Ncr ++;
     606                continue;
     607            }
     608
     609            // likely unsaturated extended source (too large to be stellar)
     610            if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) {
     611                source->type = PM_SOURCE_TYPE_EXTENDED;
     612                Next ++;
     613                continue;
     614            }
     615
     616            // the rest are probable stellar objects
     617            starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
     618            starsn_moments->n ++;
     619            starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
     620            starsn_peaks->n ++;
     621            Nstar ++;
     622
     623            // PSF star (within 1.5 sigma of clump center, S/N > limit)
     624            psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
     625            if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
     626                source->type = PM_SOURCE_TYPE_STAR;
     627                source->mode |= PM_SOURCE_MODE_PSFSTAR;
     628                Npsf ++;
     629                continue;
     630            }
    624631        }
    625632
    … …  
    636643
    637644        if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) {
    638             // Don't care about this error
    639             psErrorClear();
     645            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     646            psFree (stats);
     647            psFree (starsn_peaks);
     648            return false;
    640649        }
    641650        psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats->min, stats->max);
    … …  
    648657        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    649658        if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) {
    650             // Don't care about this error
    651             psErrorClear();
     659            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     660            psFree (stats);
     661            psFree (starsn_peaks);
     662            return false;
    652663        }
    653664        psLogMsg ("psModules.objects", 3, "SN range (peaks)  : %f - %f (%ld)\n",

  • branches/pap/psModules/src/objects/pmSource.h

    r23487 r25027  
    213213 */
    214214bool pmSourceMoments(
    215     pmSource *source,   ///< The input pmSource for which moments will be computed
    216     float radius   ///< Use a circle of pixels around the peak
     215    pmSource *source, ///< The input pmSource for which moments will be computed
     216    float radius,     ///< Use a circle of pixels around the peak
     217    float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
     218    float minSN       ///< minimum pixel significance
    217219);
    218220

  • branches/pap/psModules/src/objects/pmSourceContour.c

    r23187 r25027  
    273273            x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1);
    274274            if (x0 == x1) {
    275                 // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1);
     275              // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1);
    276276                goto pt1;
    277277            }
    … …  
    282282            x1 = findContourNeg (image, threshold, xR, yR, RIGHT);
    283283        }
    284         // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1);
     284        // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1);
    285285
    286286        xVec->data.F32[Npt + 0] = image->col0 + x0;
    … …  
    288288        yVec->data.F32[Npt + 0] = image->row0 + yR;
    289289        yVec->data.F32[Npt + 1] = image->row0 + yR;
     290
    290291        Npt += 2;
    291292
    … …  
    307308            x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1);
    308309            if (x0 == x1) {
    309                 // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1);
     310              // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1);
    310311                goto pt2;
    311312            }
    … …  
    322323        yVec->data.F32[Npt + 0] = image->row0 + yR;
    323324        yVec->data.F32[Npt + 1] = image->row0 + yR;
     325
    324326        Npt += 2;
    325327
    … …  
    334336    yVec->n = Npt;
    335337
    336     // fprintf (stderr, "done\n");
    337338    psArray *tmpArray = psArrayAlloc(2);
    338339

  • branches/pap/psModules/src/objects/pmSourceFitModel.c

    r23187 r25027  
    172172    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc);
    173173    for (int i = 0; i < dparams->n; i++) {
    174         if (psTraceGetLevel("psModules.objects") >= 4) {
    175             fprintf (stderr, "%f ", params->data.F32[i]);
    176         }
    177174        if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
    178175            continue;
    179176        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     177        if (psTraceGetLevel("psModules.objects") >= 4) {
     178            fprintf (stderr, "%f +/- %f\n", params->data.F32[i], dparams->data.F32[i]);
     179        }
    180180    }
    181181    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);

  • branches/pap/psModules/src/objects/pmSourceIO.c

    r22699 r25027  
    496496                status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
    497497            }
     498            if (!strcmp (exttype, "PS1_V2")) {
     499                status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
     500            }
    498501            if (xsrcname) {
    499502              if (!strcmp (exttype, "PS1_DEV_1")) {
    … …  
    506509                  status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
    507510              }
     511              if (!strcmp (exttype, "PS1_V2")) {
     512                  status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
     513              }
    508514            }
    509515            if (xfitname) {
    … …  
    517523                  status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
    518524              }
     525              if (!strcmp (exttype, "PS1_V2")) {
     526                  status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
     527              }
    519528            }
    520529            if (!status) {
    … …  
    562571
    563572    // not needed if only one chip
    564     if (file->fpa->chips->n == 1) return true;
    565 
     573    if (file->fpa->chips->n == 1) {
     574        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
     575        return true;
     576    }
    566577
    567578    // find the FPA phu
    … …  
    667678    psFree (outhead);
    668679
     680    pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
    669681    return true;
    670682}
    … …  
    679691
    680692    pmFPA *fpa = file->fpa;
     693
     694    pmSourceIO_ReadMatchedRefs (file->fits, fpa, config);
    681695
    682696    if (view->chip == -1) {
    … …  
    941955                sources = pmSourcesRead_CMF_PS1_V1 (file->fits, hdu->header);
    942956            }
     957            if (!strcmp (exttype, "PS1_V2")) {
     958                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
     959            }
    943960        }
    944961
    … …  
    10531070}
    10541071
     1072   

  • branches/pap/psModules/src/objects/pmSourceIO.h

    r21516 r25027  
    3939bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname);
    4040
     41bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
     42bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
     43bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
     44
    4145bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    4246
    … …  
    4852psArray *pmSourcesRead_PS1_CAL_0 (psFits *fits, psMetadata *header);
    4953psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
     54psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
    5055
    5156bool pmSourcesWritePSFs (psArray *sources, char *filename);
    … …  
    7580bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos);
    7681
     82bool pmSourceIO_WriteMatchedRefs (psFits *fits, pmFPA *fpa, pmConfig *config);
     83bool pmSourceIO_ReadMatchedRefs (psFits *fits, pmFPA *fpa, const pmConfig *config);
     84
    7785/// @}
    7886# endif /* PM_SOURCE_IO_H */

  • branches/pap/psModules/src/objects/pmSourceIO_RAW.c

    r20937 r25027  
    254254    }
    255255
     256    fprintf (f, "# %5s %5s  %8s  %7s %7s  %6s %6s  %10s %7s %7s %7s  %4s %4s %5s\n",
     257             "x", "y", "peak", "Mx", "My", "Mxx", "Myy", "Sum", "Peak", "Sky", "SN", "nPix", "type", "mode");
     258
    256259    for (i = 0; i < sources->n; i++) {
    257260        source = sources->data[i];
    258261        if (source->moments == NULL)
    259262            continue;
    260         fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
     263        fprintf (f, "%5d %5d  %8.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
    261264                 source->peak->x, source->peak->y, source->peak->value,
    262265                 source->moments->Mx, source->moments->My,

  • branches/pap/psModules/src/objects/pmSourceMasks.c

    r23184 r25027  
    4444    ADD_MASK(header, RADIAL_FLUX     , "Calculated radial flux measurements");
    4545    ADD_MASK(header, SIZE_SKIPPED    , "Could not be determine size");
     46    ADD_MASK(header, ON_SPIKE        , "Source lands in bright star spike");
     47    ADD_MASK(header, ON_GHOST        , "Source lands in bright star ghost / glint");
     48    ADD_MASK(header, OFF_CHIP        , "Source centroid lands off chip edge");
    4649
    4750    return true;

  • branches/pap/psModules/src/objects/pmSourceMasks.h

    r23184 r25027  
    3434    PM_SOURCE_MODE_RADIAL_FLUX      = 0x08000000, ///< radial flux measurements calculated
    3535    PM_SOURCE_MODE_SIZE_SKIPPED     = 0x10000000, ///< size could not be determined
     36    PM_SOURCE_MODE_ON_SPIKE         = 0x20000000, ///< peak lands on diffraction spike
     37    PM_SOURCE_MODE_ON_GHOST         = 0x40000000, ///< peak lands on ghost or glint
     38    PM_SOURCE_MODE_OFF_CHIP         = 0x80000000, ///< peak lands off edge of chip
    3639} pmSourceMode;
    3740

  • branches/pap/psModules/src/objects/pmSourceMatch.c

    r23241 r25027  
    4545                         psVector **x, psVector **y, // Coordinate vectors to return
    4646                         psVector **mag, psVector **magErr, // Magnitude and error vectors to return
     47                         psVector **indices, // Indices for sources
    4748                         const psArray *sources // Input sources
    4849    )
    … …  
    5859    *mag = psVectorRecycle(*mag, numSources, PS_TYPE_F32);
    5960    *magErr = psVectorRecycle(*magErr, numSources, PS_TYPE_F32);
     61    *indices = psVectorRecycle(*indices, numSources, PS_TYPE_S32);
    6062    float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources
    6163    long num = 0;                       // Number of valid sources
    6264    for (long i = 0; i < numSources; i++) {
    6365        pmSource *source = sources->data[i]; // Source of interest
    64         if (!source) continue;
    65         if (source->mode & SOURCE_MASK) continue;
    66         if (!isfinite(source->psfMag)) continue;
    67         if (!isfinite(source->errMag)) continue;
    68         if (source->psfMag > SOURCE_FAINTEST) continue;
     66        if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
     67            !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) {
     68            continue;
     69        }
    6970
    7071        float xSrc, ySrc;               // Coordinates of source
    … …  
    7879        (*y)->data.F32[num] = ySrc;
    7980        (*mag)->data.F32[num] = source->psfMag;
    80         (*magErr)->data.F32[num] = source->errMag ;
     81        (*magErr)->data.F32[num] = source->errMag;
     82        (*indices)->data.S32[num] = i;
    8183        num++;
    8284    }
    … …  
    8587    (*mag)->n = num;
    8688    (*magErr)->n = num;
     89    (*indices)->n = num;
    8790
    8891    if (*bounds) {
    … …  
    175178        psVector *xImage = NULL, *yImage = NULL; // Coordinates of sources
    176179        psVector *magImage = NULL, *magErrImage = NULL; // Magnitude and mag
    177 
    178         int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage,
     180        psVector *indices = NULL;     // Indices for sources
     181
     182        int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, &indices,
    179183                                      sources); // Number of sources
    180184
    … …  
    187191            for (int j = 0; j < numSources; j++) {
    188192                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    189                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     193                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     194                                 i, indices->data.S32[j]);
    190195                matches->data[j] = match;
    191196            }
    … …  
    206211            for (int j = 0, k = numMaster; j < numSources; j++, k++) {
    207212                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    208                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     213                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     214                                 i, indices->data.S32[j]);
    209215                matches->data[k] = match;
    210216            }
    … …  
    231237                    // Record the match
    232238                    pmSourceMatch *match = matches->data[index]; // Match data
    233                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     239                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     240                                     i, indices->data.S32[j]);
    234241                    numMatch++;
    235242                } else {
    236243                    // Add to the master list
    237244                    pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    238                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     245                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     246                                     i, indices->data.S32[j]);
    239247                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
    240248                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
    … …  
    469477
    470478    if (!psVectorStats(stats, trans, NULL, badImage, 0xFF)) {
     479        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
     480        psFree(stats);
     481        return -1;
     482    }
     483    // XXX handle this case better:
     484    if (isnan(stats->clippedMean))  {
    471485        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
    472486        psFree(stats);

  • branches/pap/psModules/src/objects/pmSourceMoments.c

    r21363 r25027  
    5757# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    5858
    59 bool pmSourceMoments(pmSource *source, psF32 radius)
     59bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN)
    6060{
    6161    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    8484    psF32 Y1 = 0.0;
    8585    psF32 R2 = PS_SQR(radius);
     86    psF32 minSN2 = PS_SQR(minSN);
     87    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
    8688
    8789    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    9799
    98100    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     101
     102        psF32 yDiff = row - yPeak;
     103        if (fabs(yDiff) > radius) continue;
    99104
    100105        psF32 *vPix = source->pixels->data.F32[row];
    … …  
    113118
    114119            psF32 xDiff = col - xPeak;
    115             psF32 yDiff = row - yPeak;
     120            if (fabs(xDiff) > radius) continue;
    116121
    117122            // radius is just a function of (xDiff, yDiff)
    … …  
    121126            psF32 wDiff = *vWgt;
    122127
    123             // XXX EAM : should this limit be user-defined?
    124             if (PS_SQR(pDiff) < wDiff) continue;
    125 
    126             Var += wDiff;
    127             Sum += pDiff;
    128 
    129             psF32 xWght = xDiff * pDiff;
    130             psF32 yWght = yDiff * pDiff;
    131 
    132             X1  += xWght;
    133             Y1  += yWght;
     128            // skip pixels below specified significance level
     129            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     130            if (pDiff < 0) continue;
     131
     132            if (sigma > 0.0) {
     133              // apply a pseudo-gaussian weight
     134
     135              // XXX a lot of extra flops; can we do pre-calculate?
     136              psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     137              assert (z >= 0.0);
     138              psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
     139              psF32 weight  = 1.0 / t;
     140
     141              // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
     142
     143              wDiff *= weight;
     144              pDiff *= weight;
     145            }
     146
     147            Var += wDiff;
     148            Sum += pDiff;
     149
     150            psF32 xWght = xDiff * pDiff;
     151            psF32 yWght = yDiff * pDiff;
     152
     153            X1  += xWght;
     154            Y1  += yWght;
     155
     156            // fprintf (stderr, " : %6.1f %6.1f  %8.1f %8.1f\n", X1, Y1, Sum, Var);
    134157
    135158            peakPixel = PS_MAX (*vPix, peakPixel);
    … …  
    188211    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    189212
     213        psF32 yDiff = row - yCM;
     214        if (fabs(yDiff) > radius) continue;
     215
    190216        psF32 *vPix = source->pixels->data.F32[row];
    191217        psF32 *vWgt = source->variance->data.F32[row];
    … …  
    203229
    204230            psF32 xDiff = col - xCM;
    205             psF32 yDiff = row - yCM;
     231            if (fabs(xDiff) > radius) continue;
    206232
    207233            // radius is just a function of (xDiff, yDiff)
    … …  
    215241            // XXX EAM : should this limit be user-defined?
    216242
    217             if (PS_SQR(pDiff) < wDiff) continue;
     243            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    218244            if (pDiff < 0) continue;
     245
     246            if (sigma > 0.0) {
     247              // apply a pseudo-gaussian weight
     248
     249              // XXX a lot of extra flops; can we do pre-calculate?
     250              psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     251              assert (z >= 0.0);
     252              psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
     253              psF32 weight  = 1.0 / t;
     254
     255              // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
     256
     257              wDiff *= weight;
     258              pDiff *= weight;
     259            }
    219260
    220261            Sum += pDiff;

  • branches/pap/psModules/src/objects/pmSourceUtils.c

    r23187 r25027  
    9494    source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE);
    9595
    96     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    97     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    98     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    99     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
     96    float xReadout, yReadout;
    10097
    101     // XXX fix the binning : currently not selected from concepts
    102     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    103     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    104     int xBin = 1;
    105     int yBin = 1;
     98    // if we have information about the chip & cell, adjust the coordinates chip->cell->readout
     99    // otherwise, assume 0,0 offset and 1,1 parity
     100    if (cell) {
     101      int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
     102      int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
     103      int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
     104      int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    106105
    107     // Position on the cell
    108     float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
    109     float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
     106      // XXX fix the binning : currently not selected from concepts
     107      // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
     108      // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
     109      int xBin = 1;
     110      int yBin = 1;
    110111
    111     // Position on the readout
    112     // float xReadout = CELL_TO_READOUT(xCell, x0Readout);
    113     // float yReadout = CELL_TO_READOUT(yCell, y0Readout);
     112      // Position on the cell
     113      float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
     114      float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
     115
     116      // Position on the readout
     117      // float xReadout = CELL_TO_READOUT(xCell, x0Readout);
     118      // float yReadout = CELL_TO_READOUT(yCell, y0Readout);
     119      xReadout = xCell;
     120      yReadout = yCell;
     121    } else {
     122      xReadout = xChip;
     123      yReadout = yChip;
     124    }
    114125   
    115     pmSourceDefinePixels (source, readout, xCell, yCell, radius);
     126    pmSourceDefinePixels (source, readout, xReadout, yReadout, radius);
    116127
    117128    return (source);

  • branches/pap/psModules/src/psmodules.h

    r23751 r25027  
    4747#include <pmFPAfileDefine.h>
    4848#include <pmFPAfileFitsIO.h>
     49#include <pmFPAfileFringeIO.h>
    4950#include <pmFPAfileIO.h>
    5051#include <pmFPARead.h>
    … …  
    7677#include <pmDark.h>
    7778#include <pmRemnance.h>
     79#include <pmPattern.h>
    7880
    7981// the following headers are from psModule:astrom

  • branches/pap/psModules/test/camera/tap_pmFPAMaskW.c

    r21474 r25027  
    356356    // ----------------------------------------------------------------------
    357357    // pmReadoutSetVariance() tests: NULL inputs
    358     // bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
     358    // bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson)
    359359    if (1) {
    360360        psMemId id = psMemGetId();
    … …  
    367367
    368368        // Set readout == NULL, ensure pmReadoutSetVariance() returnes FALSE, with no seg faults, memory leaks
    369         rc = pmReadoutSetVariance(NULL, false);
    370         ok(!rc, "pmReadoutSetVariance(NULL, false) returned FALSE with null pmReadout input");
     369        rc = pmReadoutSetVariance(NULL, NULL, false);
     370        ok(!rc, "pmReadoutSetVariance(NULL, NULL, false) returned FALSE with null pmReadout input");
    371371
    372372
    … …  
    375375        rc|= psMetadataAddF32(readout->parent->concepts, PS_LIST_HEAD, "CELL.READNOISE", PS_META_REPLACE, NULL, CELL_READNOISE);
    376376        ok(rc, "Set GAIN and READNOISE in cell->concepts successfully");
    377 
     377/*
     378 * Getting the section below to run requires generating a noiseMap
     379 *
    378380        // Call pmReadoutSetVariance() and then verify that the mask data was set correctly
    379381        rc = pmReadoutSetVariance(readout, false);
    … …  
    419421
    420422        ok(!errorFlag, "pmReadoutSetWeight() set the weight values correctly (Poisson)");
     423*/     
    421424        psFree(fpa);   
    422425        psFree(camera);

  • branches/pap/psModules/test/concepts/tap_pmConcepts.c

    r21223 r25027  
    7777}
    7878
     79psMetadataItem *dummyConceptCopier(
     80    const psMetadataItem *source,
     81    const psMetadataItem *target,
     82    const psMetadata *cameraFormat,
     83    const pmFPA *fpa,
     84    const pmChip *chip,
     85    const pmCell *cell)
     86{
     87    if (target == NULL ||
     88        source == PM_CONCEPT_SOURCE_NONE ||
     89        cameraFormat == NULL ||
     90        fpa == NULL ||
     91        chip == NULL ||
     92        cell == NULL) {
     93        printf("dummyConceptCopier() args are NULL\n");
     94    }
     95    return(NULL);
     96}
    7997
    8098int main(int argc, char* argv[])
    … …  
    92110        psMetadataItem *blank = psMetadataItemAlloc("myItem1", PS_DATA_BOOL, "I am a boolean", true);
    93111        pmConceptSpec *tmp = pmConceptSpecAlloc(blank, dummyConceptParser,
    94             dummyConceptFormatter, true);
     112            dummyConceptFormatter, dummyConceptCopier, true);
    95113        ok(tmp != NULL, "pmConceptSpecAlloc() returned non-NULL");
    96114        skip_start(tmp == NULL, 4, "Skipping tests because pmConceptSpecAlloc() returned NULL");
    … …  
    98116        ok(tmp->parse == dummyConceptParser, "pmConceptSpecAlloc() set the ->parse member correctly");
    99117        ok(tmp->format == dummyConceptFormatter, "pmConceptSpecAlloc() set the ->format member correctly");
     118        ok(tmp->copy == dummyConceptCopier, "pmConceptSpecAlloc() set the ->copy member correctly");
    100119        ok(tmp->required == true, "pmConceptSpecAlloc() set the ->required member correctly");
    101120        skip_end();
    … …  
    109128    {
    110129        psMemId id = psMemGetId();
    111         pmConceptSpec *tmp = pmConceptSpecAlloc(NULL, NULL, NULL, false);
     130        pmConceptSpec *tmp = pmConceptSpecAlloc(NULL, NULL, NULL, NULL, false);
    112131        ok(tmp != NULL, "pmConceptSpecAlloc() returned non-NULL with NULL inputs");
    113132        psFree(tmp);

  • branches/pap/psModules/test/objects/tap_pmGrowthCurve.c

    r21223 r25027  
    129129
    130130        source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1;
    131         source->mode = PM_SOURCE_MODE_SUBTRACTED;
     131        source->mode = PM_SOURCE_MODE_PSFSTAR;
    132132
    133133        source->modelPSF->radiusFit = 15.0;
    … …  
    232232
    233233        source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1;
    234         source->mode = PM_SOURCE_MODE_SUBTRACTED;
     234        source->mode = PM_SOURCE_MODE_PSFSTAR;
    235235
    236236        source->modelPSF->radiusFit = 15.0;
    … …  
    334334
    335335        source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1;
    336         source->mode = PM_SOURCE_MODE_SUBTRACTED;
     336        source->mode = PM_SOURCE_MODE_PSFSTAR;
    337337
    338338        source->modelPSF->radiusFit = 15.0;

  • branches/pap/psModules/test/objects/tap_pmSource.c

    r21471 r25027  
    432432        psArray *sources = psArrayAlloc(NUM_SOURCES);
    433433        psMetadata *recipe = psMetadataAlloc();
    434         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     434        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    435435        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    436436        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);
    … …  
    452452        psArray *sources = psArrayAlloc(NUM_SOURCES);
    453453        psMetadata *recipe = psMetadataAlloc();
    454         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     454        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    455455        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    456456        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);
    … …  
    490490        }
    491491        psMetadata *recipe = psMetadataAlloc();
    492         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     492        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    493493        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    494494        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);

  • branches/pap/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c

    r21223 r25027  
    4848        psMetadata *tableHeader = psMetadataAlloc();
    4949        psString extname = psStringCopy("ext");
    50         bool rc = pmSourcesWrite_PS1_DEV_1(NULL, sources, imageHeader, tableHeader, extname, NULL);
     50        bool rc = pmSourcesWrite_PS1_DEV_1(NULL, sources, imageHeader, tableHeader, extname);
    5151        ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL psFits input parameter");
    5252        psFree(fitsFile);
    … …  
    7272        psMetadata *tableHeader = psMetadataAlloc();
    7373        psString extname = psStringCopy("ext");
    74         bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, NULL, imageHeader, tableHeader, extname, NULL);
     74        bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, NULL, imageHeader, tableHeader, extname);
    7575        ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL pmSource input parameter");
    7676        psFree(fitsFile);
    … …  
    9696        psMetadata *tableHeader = psMetadataAlloc();
    9797        psString extname = psStringCopy("ext");
    98         bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, NULL, NULL);
     98        bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, NULL);
    9999        ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL extname input parameter");
    100100        psFree(fitsFile);
    … …  
    216216        psMetadata *tableHeader = psMetadataAlloc();
    217217        psString extname = psStringCopy("ext");
    218         bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, extname, NULL);
     218        bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, extname);
    219219        ok(rc == true, "pmSourcesWrite_PS1_DEV_1() returned TRUE with acceptable input parameters");
    220220        psFree(fitsFile);
Note: See TracChangeset for help on using the changeset viewer.