IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 26, 2009, 1:59:32 PM (17 years ago)
Author:
beaumont
Message:

merged with trunk

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
63 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/psModules

  • branches/cnb_branches/cnb_branch_20090301/psModules/configure.ac

    r14187 r24244  
    305305echo "PSMODULES_CFLAGS: $PSMODULES_CFLAGS"
    306306echo "PSMODULE_LIBS: $PSMODULES_LIBS"
     307
     308IPP_VERSION(PSMODULES)
     309AM_CONDITIONAL([HAVE_SVNVERSION], [test "x$SVNVERSION" != x])
     310AM_CONDITIONAL([HAVE_SVN], [test "x$SVN" != x])
    307311
    308312dnl ------- enable -Werror after all of the probes have run ---------

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryDistortion.c

    r23594 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryObjects.c

    r21422 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryObjects.h

    r20801 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryUtils.c

    r19991 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/astrom/pmAstrometryWCS.c

    r21430 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPA.c

    r21363 r24244  
    298298
    299299
    300 pmCell *pmCellAlloc(
    301     pmChip *chip,
    302     const char *name)
     300pmCell *pmCellAlloc(pmChip *chip, const char *name)
    303301{
    304302    pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell));
    … …  
    323321    }
    324322    tmpCell->conceptsRead = PM_CONCEPT_SOURCE_NONE;
    325     // XXX does this work?  moved to conceptsRead... pmConceptsBlankCell(tmpCell);
     323    pmConceptsBlankCell(tmpCell);
    326324
    327325    return tmpCell;
    … …  
    362360    }
    363361    tmpChip->conceptsRead = PM_CONCEPT_SOURCE_NONE;
    364     // XXX does this work?  moved to conceptsRead... pmConceptsBlankChip(tmpChip);
     362    pmConceptsBlankChip(tmpChip);
    365363    return tmpChip;
    366364}
    … …  
    393391    }
    394392    tmpFPA->conceptsRead = PM_CONCEPT_SOURCE_NONE;
    395     // XXX does this work?  moved to conceptsRead... pmConceptsBlankFPA(tmpFPA);
     393    pmConceptsBlankFPA(tmpFPA);
    396394
    397395    // this may be somewhat pedantic, but it makes these things consistent

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPABin.c

    r21183 r24244  
    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
    … …  
    5865            }
    5966
    60             outImage->data.F32[yOut][xOut] = numPix > 0 ? sum / numPix : NAN;
     67            float imageValue, maskValue;// Values to set
     68            if (numPix > 0) {
     69                imageValue = sum / numPix;
     70                maskValue = 0;
     71            } else {
     72                imageValue = NAN;
     73                maskValue = maskVal;
     74            }
     75            outImage->data.F32[yOut][xOut] = imageValue;
     76            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
    6177            xStart = xStop;
    6278        }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPACopy.c

    r23351 r24244  
    399399    if (targetHDU) {
    400400        pmHDU *sourceHDU = pmHDUFromCell(source); // The source HDU
    401         if (!sourceHDU) {
    402             psWarning("Unable to copy header: no source header.");
    403         } else if (sourceHDU->header) {
     401        if (sourceHDU && sourceHDU->header) {
    404402            targetHDU->header = psMetadataCopy(targetHDU->header, sourceHDU->header);
    405403        }
    … …  
    481479    if (targetFPA && sourceFPA) {
    482480        if (targetFPA->toSky) {
    483             psAssert (targetFPA->toSky == sourceFPA->toSky, "chips within FPA have inconsistent astrometry references");
    484         } else {
    485             targetFPA->toSky = psMemIncrRefCounter (sourceFPA->toSky);
    486         }
     481            psFree(targetFPA->toSky);
     482        }
     483        targetFPA->toSky = psMemIncrRefCounter (sourceFPA->toSky);
    487484        if (targetFPA->toTPA) {
    488             psAssert (targetFPA->toTPA == sourceFPA->toTPA, "chips within FPA have inconsistent astrometry references");
    489         } else {
    490             targetFPA->toTPA = psMemIncrRefCounter (sourceFPA->toTPA);
    491         }
     485            psFree(targetFPA->toTPA);
     486        }
     487        targetFPA->toTPA = psMemIncrRefCounter (sourceFPA->toTPA);
    492488        if (targetFPA->fromTPA) {
    493             psAssert (targetFPA->fromTPA == sourceFPA->fromTPA, "chips within FPA have inconsistent astrometry references");
    494         } else {
    495             targetFPA->fromTPA = psMemIncrRefCounter (sourceFPA->fromTPA);
    496         }
     489            psFree(targetFPA->fromTPA);
     490        }
     491        targetFPA->fromTPA = psMemIncrRefCounter (sourceFPA->fromTPA);
    497492    }
    498493

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPAMaskWeight.c

    r23351 r24244  
    482482        return false;
    483483    }
    484     float stdev = psStatsGetValue(stdevStats, stdevStat); // Stadard deviation of fluxes
     484    float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes
    485485    psFree(stdevStats);
    486486    psFree(noise);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPARead.c

    r23594 r24244  
    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
    … …  
    253251    int readdir = psMetadataLookupS32(&mdok, cell->concepts, "CELL.READDIR"); // Read direction
    254252    if (!mdok || readdir == 0 || (readdir != 1 && readdir != 2)) {
    255         psError(PS_ERR_IO, true, "CELL.READDIR is not set to -1 or +1.\n");
     253        psError(PS_ERR_IO, true, "CELL.READDIR is not set to 1 or 2.\n");
    256254        return false;
    257255    }
    … …  
    592590    }
    593591
    594     // Determine the number of scans to read
    595     int lastScan = readoutSetLastScan(readout, type, thisScan + numScans);
     592    int origThisScan = thisScan;        // Original value of thisScan (starting point for read)
    596593    if (thisScan == 0) {
    597594        overlap = 0;
    … …  
    601598        thisScan = 0;
    602599    }
    603     readoutSetThisScan(readout, type, thisScan);
    604600
    605601    // Calculate limits, adjust readout->row0,col0
    … …  
    620616        }
    621617    }
     618    int lastScan = origThisScan + numScans; // Last scan to read
     619
     620    readoutSetThisScan(readout, type, thisScan);
     621    readoutSetLastScan(readout, type, lastScan);
    622622
    623623    // Blow away existing data.

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmFPAfile.c

    r23594 r24244  
    205205    }
    206206    if (strstr(newName, "{CHIP.NAME}")) {
     207        const char *name = NULL;        // Name of chip
    207208        pmChip *chip = pmFPAviewThisChip(view, fpa);
    208209        if (chip) {
    209             char *name = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    210             if (name) {
    211                 psStringSubstitute(&newName, name, "{CHIP.NAME}");
    212             }
    213         }
     210            name = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
     211            psAssert(name, "All chips should have a name");
     212        } else {
     213            name = "fpa";
     214        }
     215        psStringSubstitute(&newName, name, "{CHIP.NAME}");
     216
    214217    }
    215218    if (strstr(newName, "{CHIP.ID}")) {
    … …  
    255258    }
    256259    if (strstr(newName, "{CELL.NAME}")) {
     260        const char *name = NULL;        // Name of cell
    257261        pmCell *cell = pmFPAviewThisCell(view, fpa);
    258262        if (cell) {
    259             char *name = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
    260             if (name) {
    261                 psStringSubstitute(&newName, name, "{CELL.NAME}");
    262             }
    263         }
     263            name = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
     264            psAssert(name, "All cells should have a name");
     265        } else {
     266            name = "chip";
     267        }
     268        psStringSubstitute(&newName, name, "{CELL.NAME}");
    264269    }
    265270    if (strstr(newName, "{CELL.N}")) {

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/camera/pmReadoutFake.c

    r21104 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/concepts/pmConcepts.c

    r22699 r24244  
    421421
    422422// Interpolate the concept.  Generalises the FPA/Chip/Cell
    423 #define CONCEPT_INTERPOLATE(SOURCE, NAME) \
     423#define CONCEPT_INTERPOLATE(SOURCE, NAME, DEFAULT) \
    424424    if (strncmp(concept, NAME, strlen(NAME)) == 0) { \
    425         if (!(SOURCE)) { \
    426             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot interpolate %s because %s not provided", \
    427                     concept, NAME); \
    428             psFree(string); \
    429             return NULL; \
    430         } \
    431         psMetadataItem *item = psMetadataLookup((SOURCE)->concepts, concept); /* Item with concept */ \
    432         if (!item) { \
    433             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Can't find concept %s in %s", concept, NAME); \
    434             psFree(string); \
    435             return NULL; \
    436         } \
    437         \
    438         psString value = psMetadataItemParseString(item); /* Value of concept */ \
    439         if (!value) { \
    440             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to parse concept %s", concept); \
    441             psFree(string); \
    442             return NULL; \
     425        psString value = NULL; /* Value of concept */ \
     426        if (SOURCE) { \
     427            psMetadataItem *item = psMetadataLookup((SOURCE)->concepts, concept); /* Item with concept */ \
     428            if (!item) { \
     429                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Can't find concept %s in %s", concept, NAME); \
     430                psFree(string); \
     431                return NULL; \
     432            } \
     433            \
     434            value = psMetadataItemParseString(item); /* Value of concept */ \
     435            if (!value) { \
     436                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to parse concept %s", concept); \
     437                psFree(string); \
     438                return NULL; \
     439            } \
     440        } else { \
     441            if (!(DEFAULT)) { \
     442                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to determine concept %s", concept); \
     443                psFree(string); \
     444                return NULL; \
     445            } \
     446            value = psStringCopy(DEFAULT); \
    443447        } \
    444448        \
    … …  
    483487        psTrace("psModules.concepts", 7, "Interpolating concept %s", concept);
    484488
    485         CONCEPT_INTERPOLATE(fpa,  "FPA");
    486         CONCEPT_INTERPOLATE(chip, "CHIP");
    487         CONCEPT_INTERPOLATE(cell, "CELL");
     489        CONCEPT_INTERPOLATE(fpa,  "FPA", NULL);
     490        CONCEPT_INTERPOLATE(chip, "CHIP", "fpa");
     491        CONCEPT_INTERPOLATE(cell, "CELL", "chip");
    488492
    489493        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised concept: %s", concept);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/concepts/pmConceptsStandard.c

    r23594 r24244  
    3333break;
    3434
     35
     36// Format type for time
     37typedef 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
     43} conceptTimeFormatType;
     44
     45// Format for time
     46typedef struct {
     47    conceptTimeFormatType format;       // Format type for time
     48    bool separate;                      // Date and time stored separately?
     49    bool pre2000;                       // Year is pre-2000 (two digits only)?
     50} conceptTimeFormat;
    3551
    3652
    … …  
    630646}
    631647
     648static conceptTimeFormat conceptGetTimeFormat(const char *name, // Concept name ("CELL.TIME" or "FPA.TIME")
     649                                              const psMetadata *cameraFormat // Camera format
     650    )
     651{
     652    conceptTimeFormat time;               // Time format, to return
     653    time.format = TIME_FORMAT_YYYYMMDD;
     654    time.separate = false;
     655    time.pre2000 = false;
     656
     657    bool mdok = true;                   // Status of MD lookup
     658    psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS"); // The formats
     659    if (mdok && formats) {
     660        psString format = psMetadataLookupStr(&mdok, formats, name); // The formats for eg, CELL.TIME
     661        if (mdok && format && strlen(format) > 0) {
     662            psList *formatList = psStringSplit(format, " ,;", false); // List of formats specified
     663            psListIterator *formatListIter = psListIteratorAlloc(formatList, PS_LIST_HEAD, false); // Iterator
     664            while ((format = psListGetAndIncrement(formatListIter))) {
     665                if (strcasecmp(format, "SEPARATE") == 0) {
     666                    time.separate = true;
     667                } else if (strcasecmp(format, "YYYYMMDD") == 0) {
     668                    time.format = TIME_FORMAT_YYYYMMDD;
     669                } else if (strcasecmp(format, "MMDDYYYY") == 0) {
     670                    time.format = TIME_FORMAT_MMDDYYYY;
     671                } else if (strcasecmp(format, "DDMMYYYY") == 0) {
     672                    time.format = TIME_FORMAT_DDMMYYYY;
     673                } else if (strcasecmp(format, "ISO") == 0) {
     674                    time.format = TIME_FORMAT_YYYYMMDD;
     675                } else if (strcasecmp(format, "YEAR.FIRST") == 0) {
     676                    time.format = TIME_FORMAT_YYYYMMDD;
     677                } else if (strcasecmp(format, "USA") == 0) {
     678                    time.format = TIME_FORMAT_MMDDYYYY;
     679                } else if (strcasecmp(format, "BACKWARDS") == 0) {
     680                    time.format = TIME_FORMAT_DDMMYYYY;
     681                } else if (strcasecmp(format, "PRE2000") == 0) {
     682                    time.pre2000 = true;
     683                } else if (strcasecmp(format, "MJD") == 0) {
     684                    time.format = TIME_FORMAT_MJD;
     685                } else if (strcasecmp(format, "JD") == 0) {
     686                    time.format = TIME_FORMAT_JD;
     687                } else {
     688                    psWarning("Unrecognised FORMATS option for %s: %s --- ignored.", name, format);
     689                }
     690            }
     691            psFree(formatListIter);
     692            psFree(formatList);
     693        }
     694    }
     695
     696    return time;
     697}
     698
    632699// Determine the corresponding TIMESYS for one of the TIME concepts
    633700static psTimeType conceptGetTimesysForTime(const char *name, // Concept name ("CELL.TIME" or "FPA.TIME")
    … …  
    734801    psTimeType timeSys = conceptGetTimesysForTime(pattern->name, fpa, chip, cell); // Time system
    735802
    736     // Work out how the time is represented
    737     bool separateTime = false;
    738     bool usaTime = false;               // Is the time specified in USA (MM-DD-YYYY) order?
    739     bool backwardsTime = false;         // Is the time specified in backwards (DD-MM-YYYY) order?
    740     bool yearFirst = false;            // Is the time specified in yearFirst (YYYY-MM-DD) order?
    741     bool pre2000Time = false;           // Is the time specified pre-2000 (where the year only has two digits)
    742     bool mjdTime = false;               // Is the time specified a MJD?
    743     bool jdTime = false;                // Is the time specified a JD?
    744 
    745     // Get format
    746     bool mdok;                          // Status of MD lookup
    747     psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS");
    748     if (!mdok || !formats) {
    749         psError(PS_ERR_UNKNOWN, true, "Unable to find FORMATS in camera configuration.\n");
    750         return NULL;
    751     }
    752 
    753     psString timeFormat = psMetadataLookupStr(&mdok, formats, pattern->name);
    754     if (!mdok || strlen(timeFormat) == 0) {
    755         psError(PS_ERR_UNKNOWN, true, "Unable to find %s in FORMATS.\n", pattern->name);
    756         return NULL;
    757     }
    758 
    759     // Parse the time format
    760     // why should more than one be allowed??
    761     psList *timeFormats = psStringSplit(timeFormat, " ,;", false); // List of the format options
    762     psListIterator *timeFormatsIter = psListIteratorAlloc(timeFormats, PS_LIST_HEAD, false); // Iter
    763     while ((timeFormat = psListGetAndIncrement(timeFormatsIter))) {
    764         if (strcasecmp(timeFormat, "SEPARATE") == 0) {
    765             separateTime = true;
    766         } else if (strcasecmp(timeFormat, "USA") == 0) {
    767             usaTime = true;
    768             backwardsTime = false;
    769             yearFirst = false;
    770             jdTime = false;
    771             mjdTime = false;
    772         } else if (strcasecmp(timeFormat, "BACKWARDS") == 0) {
    773             backwardsTime = true;
    774             usaTime = false;
    775             yearFirst = false;
    776             jdTime = false;
    777             mjdTime = false;
    778         } else if (strcasecmp(timeFormat, "YEAR.FIRST") == 0) {
    779             yearFirst = true;
    780             usaTime = false;
    781             backwardsTime = false;
    782             jdTime = false;
    783             mjdTime = false;
    784         } else if (strcasecmp(timeFormat, "PRE2000") == 0) {
    785             pre2000Time = true;
    786         } else if (strcasecmp(timeFormat, "MJD") == 0) {
    787             mjdTime = true;
    788             jdTime = false;
    789             yearFirst = false;
    790             backwardsTime = false;
    791             usaTime = false;
    792         } else if (strcasecmp(timeFormat, "JD") == 0) {
    793             jdTime = true;
    794             mjdTime = false;
    795             yearFirst = false;
    796             backwardsTime = false;
    797             usaTime = false;
    798         } else {
    799             psError(PS_ERR_UNKNOWN, true, "Unrecognised FORMATS option for %s: %s --- "
    800                     "ignored.\n", pattern->name, timeFormat);
    801             psFree(timeFormatsIter);
    802             psFree(timeFormats);
    803             return NULL;
    804         }
    805     }
    806     psFree(timeFormatsIter);
    807     psFree(timeFormats);
     803    conceptTimeFormat timeFormat = conceptGetTimeFormat(pattern->name, cameraFormat); // Format for time
    808804
    809805    psTime *time = NULL;                // The time
    810806    switch (concept->type) {
    811807      case PS_DATA_LIST: {
    812           if (!separateTime) {
     808          if (!timeFormat.separate) {
    813809              psWarning ("DATE and TIME stored separately, but not specified in format\n");
    814810          }
    … …  
    837833              return NULL;
    838834          }
    839           if (backwardsTime) {
    840               // Need to switch days and years
    841               int temp = day;
    842               day = year;
    843               year = temp;
    844           }
    845           if (usaTime) {
    846               // Need to switch everything around.... Yanks!
    847               int temp = day;
    848               day = month;
    849               month = year;
    850               year = temp;
     835          switch (timeFormat.format) {
     836            case TIME_FORMAT_DDMMYYYY: {
     837                // Need to switch days and years
     838                int temp = day;
     839                day = year;
     840                year = temp;
     841                break;
     842            }
     843            case TIME_FORMAT_MMDDYYYY: {
     844                // Need to switch everything around.... Yanks!
     845                int temp = day;
     846                day = month;
     847                month = year;
     848                year = temp;
     849                break;
     850            }
     851            default:
     852              break;
    851853          }
    852854          if (year < 100) {
    853               if (pre2000Time) {
     855              if (timeFormat.pre2000) {
    854856                  year += 1900;
    855857              } else {
    … …  
    898900      case PS_DATA_STRING: {
    899901          psString timeString = concept->data.V;   // String with the time
    900           if (jdTime) {
    901               double timeValue = strtod (timeString, NULL);
    902               time = psTimeFromJD(timeValue);
    903           } else if (mjdTime) {
    904               double timeValue = strtod (timeString, NULL);
    905               time = psTimeFromMJD(timeValue);
    906           } else {
    907               // It's ISO
    908               time = psTimeFromISO(timeString, timeSys);
    909           } // Interpreting the time string
     902          switch (timeFormat.format) {
     903            case TIME_FORMAT_JD: {
     904                double timeValue = strtod (timeString, NULL);
     905                time = psTimeFromJD(timeValue);
     906                break;
     907            }
     908            case TIME_FORMAT_MJD: {
     909                double timeValue = strtod (timeString, NULL);
     910                time = psTimeFromMJD(timeValue);
     911                break;
     912            }
     913            case TIME_FORMAT_YYYYMMDD: {
     914                // this is ISO-standard
     915                time = psTimeFromISO(timeString, timeSys);
     916                break;
     917            }
     918            default:
     919              psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to interpret time string: %s", timeString);
     920              return NULL;
     921          }
    910922          break;
    911923      }
    912924      case PS_TYPE_F32: {
    913925          double timeValue = (double)concept->data.F32;
    914           if (jdTime) {
     926          switch (timeFormat.format) {
     927            case TIME_FORMAT_JD:
    915928              time = psTimeFromJD(timeValue);
    916           } else if (mjdTime) {
     929              break;
     930            case TIME_FORMAT_MJD:
    917931              time = psTimeFromMJD(timeValue);
    918           } else {
    919               psError(PS_ERR_UNKNOWN, true, "Not sure how to parse %s (%f) --- trying JD\n",
    920                       pattern->name, timeValue);
    921               time = psTimeFromJD(timeValue);
     932              break;
     933            default:
     934              psError(PS_ERR_UNKNOWN, true, "Unable to interpret time %s (%f)", pattern->name, timeValue);
     935              return NULL;
    922936          }
    923937          break;
    … …  
    925939      case PS_TYPE_F64: {
    926940          double timeValue = (double)concept->data.F64;
    927           if (jdTime) {
     941          switch (timeFormat.format) {
     942            case TIME_FORMAT_JD:
    928943              time = psTimeFromJD(timeValue);
    929           } else if (mjdTime) {
     944              break;
     945            case TIME_FORMAT_MJD:
    930946              time = psTimeFromMJD(timeValue);
    931           } else {
    932               psError(PS_ERR_UNKNOWN, true, "Not sure how to parse %s (%f) --- trying JD\n",
    933                       pattern->name, timeValue);
    934               time = psTimeFromJD(timeValue);
     947              break;
     948            default:
     949              psError(PS_ERR_UNKNOWN, true, "Unable to interpret time %s (%f)", pattern->name, timeValue);
     950              return NULL;
    935951          }
    936952          break;
    … …  
    946962    }
    947963
    948     if (jdTime || mjdTime) {
     964    // Set the time system appropriately
     965    switch (timeFormat.format) {
     966      case TIME_FORMAT_JD:
     967      case TIME_FORMAT_MJD:
    949968        conceptSetTimesysForTime(pattern->name, fpa, chip, cell, PS_TIME_TAI);
    950     } else {
     969        break;
     970      default:
    951971        time->type = timeSys;
     972        break;
    952973    }
    953974
    … …  
    11761197                                          const pmCell *cell)
    11771198{
     1199    psString timeName = psStringCopy(concept->name); // Name of corresponding TIME concept
     1200    psStringSubstitute(&timeName, "TIME", "TIMESYS");
     1201
     1202    conceptTimeFormat timeFormat = conceptGetTimeFormat(timeName, cameraFormat); // Format for time
     1203    psFree(timeName);
     1204
     1205    psTimeType timesys = concept->data.S32; // Time system
     1206
     1207    // JD and MJD are converted to TAI before writing
     1208    switch (timeFormat.format) {
     1209      case TIME_FORMAT_JD:
     1210      case TIME_FORMAT_MJD:
     1211        timesys = PS_TIME_TAI;
     1212        break;
     1213      default:
     1214        break;
     1215    }
     1216
    11781217    psString sys = NULL;            // String to store
    1179     switch (concept->data.S32) {
     1218    switch (timesys) {
    11801219      case PS_TIME_TAI:
    11811220        sys = psStringCopy("TAI");
    … …  
    12111250    psTimeConvert(time, timeSys);
    12121251
    1213     // Work out the format
    1214     bool separateTime = false;          // Are the date and time stored separately?
    1215     bool pre2000Time = false;           // Is the year in pre-2000 format (two digits only)?
    1216     bool backwardsTime = false;         // Is the date stored backwards (DD-MM-YYYY)?
    1217     bool usaTime = false;               // Is the date stored in USA order (MM-DD-YYYY)?
    1218     bool jdTime = false;                // Is the date stored as a JD?
    1219     bool mjdTime = false;               // Is the date stored as a MJD?
    1220     bool yearFirst = false;
    1221 
    1222     bool mdok = true;                   // Status of MD lookup
    1223     psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS"); // The formats
    1224     if (mdok && formats) {
    1225         psString format = psMetadataLookupStr(&mdok, formats, concept->name); // The formats for eg, CELL.TIME
    1226         if (mdok && format && strlen(format) > 0) {
    1227             psList *formatList = psStringSplit(format, " ,;", false); // List of formats specified
    1228             psListIterator *formatListIter = psListIteratorAlloc(formatList, PS_LIST_HEAD, false); // Iterator
    1229             while ((format = psListGetAndIncrement(formatListIter))) {
    1230                 if (strcasecmp(format, "SEPARATE") == 0) {
    1231                     separateTime = true;
    1232                 } else if (strcasecmp(format, "USA") == 0) {
    1233                     usaTime = true;
    1234                     backwardsTime = false;
    1235                     jdTime = false;
    1236                     mjdTime = false;
    1237                 } else if (strcasecmp(format, "BACKWARDS") == 0) {
    1238                     backwardsTime = true;
    1239                     usaTime = false;
    1240                     mjdTime = false;
    1241                     jdTime = false;
    1242                 } else if (strcasecmp(format, "YEAR.FIRST") == 0) {
    1243                     yearFirst = true;
    1244                     usaTime = false;
    1245                     backwardsTime = false;
    1246                     jdTime = false;
    1247                     mjdTime = false;
    1248                 } else if (strcasecmp(format, "PRE2000") == 0) {
    1249                     pre2000Time = true;
    1250                 } else if (strcasecmp(format, "MJD") == 0) {
    1251                     mjdTime = true;
    1252                     usaTime = false;
    1253                     backwardsTime = false;
    1254                     jdTime = false;
    1255                     separateTime = false;
    1256                 } else if (strcasecmp(format, "JD") == 0) {
    1257                     jdTime = true;
    1258                     usaTime = false;
    1259                     backwardsTime = false;
    1260                     mjdTime = false;
    1261                     separateTime = false;
    1262                 } else {
    1263                     psWarning("Unrecognised FORMATS option for %s: %s --- "
    1264                               "ignored.\n", concept->name, format);
    1265                 }
    1266             }
    1267             psFree(formatListIter);
    1268             psFree(formatList);
    1269         }
    1270     }
    1271 
    1272     if (separateTime) {
     1252    conceptTimeFormat timeFormat = conceptGetTimeFormat(concept->name, cameraFormat); // Format for time
     1253
     1254    if (timeFormat.separate) {
    12731255        // We're working with two separate headers --- construct a list with the date and time separately
    12741256        psString dateTimeString = psTimeToISO(time); // String representation
    … …  
    12801262        // Need to format the strings....
    12811263        // XXX: Couldn't be bothered doing these right now
    1282         if (pre2000Time) {
     1264        if (timeFormat.pre2000) {
    12831265            psError(PS_ERR_UNKNOWN, true, "Don't you realise it's the twenty-first century?\n");
    12841266            return NULL;
    12851267        }
    1286         if (backwardsTime) {
    1287             int day, month, year;
    1288             psTrace ("psModules.concepts", 5, "ISO time has year first, convert to DD-MM-YYYY");
    1289             sscanf (dateString, "%d-%d-%d", &year, &month, &day);
    1290             sprintf (dateString, "%02d-%02d-%04d", day, month, year);
    1291             // XXX fix this for str length
    1292         }
    1293         if (usaTime) {
    1294             int day, month, year;
    1295             psTrace ("psModules.concepts", 5, "ISO time has year first, convert to MM-DD-YYYY");
    1296             sscanf (dateString, "%d-%d-%d", &year, &month, &day);
    1297             sprintf (dateString, "%02d-%02d-%04d", month, day, year);
    1298             // XXX fix this for str length
    1299         }
    1300         if (yearFirst) {
    1301             psTrace ("psModules.concepts", 5, "ISO time has year first, no adjustment needed");
     1268
     1269        switch (timeFormat.format) {
     1270          case TIME_FORMAT_DDMMYYYY: {
     1271              int day, month, year;
     1272              psTrace ("psModules.concepts", 5, "ISO time has year first, convert to DD-MM-YYYY");
     1273              sscanf (dateString, "%d-%d-%d", &year, &month, &day);
     1274              sprintf (dateString, "%02d-%02d-%04d", day, month, year);
     1275              // XXX fix this for str length
     1276              break;
     1277          }
     1278          case TIME_FORMAT_MMDDYYYY: {
     1279              int day, month, year;
     1280              psTrace ("psModules.concepts", 5, "ISO time has year first, convert to MM-DD-YYYY");
     1281              sscanf (dateString, "%d-%d-%d", &year, &month, &day);
     1282              sprintf (dateString, "%02d-%02d-%04d", month, day, year);
     1283              // XXX fix this for str length
     1284              break;
     1285          }
     1286          default:
     1287            break;
    13021288        }
    13031289
    … …  
    13131299        psListAdd(dateTime, PS_LIST_TAIL, timeItem);
    13141300
    1315         psMetadataItem *item = psMetadataItemAllocPtr(concept->name, PS_DATA_LIST, concept->comment,
    1316                                                       dateTime);
    1317         psFree (dateItem);
    1318         psFree (timeItem);
    1319         psFree (dateTime);
     1301        psMetadataItem *item = psMetadataItemAllocPtr(concept->name, PS_DATA_LIST,
     1302                                                      concept->comment, dateTime);
     1303        psFree(dateItem);
     1304        psFree(timeItem);
     1305        psFree(dateTime);
    13201306        return item;
    13211307    }
    1322     if (jdTime) {
    1323         double jd = psTimeToMJD(time);
    1324         return psMetadataItemAllocF64(concept->name, concept->comment, jd);
    1325     }
    1326     if (mjdTime) {
    1327         double mjd = psTimeToMJD(time);
    1328         return psMetadataItemAllocF64(concept->name, concept->comment, mjd);
    1329     }
    1330 
    1331     // If we've gotten this far, so it's straight ISO.
    1332     psString dateTimeString = psTimeToISO(time); // String representation
    1333     psMetadataItem *item = psMetadataItemAllocStr(concept->name, concept->comment, dateTimeString);
    1334     psFree(dateTimeString);
    1335     return item;
     1308
     1309    switch (timeFormat.format) {
     1310      case TIME_FORMAT_JD: {
     1311          double jd = psTimeToMJD(time);
     1312          return psMetadataItemAllocF64(concept->name, concept->comment, jd);
     1313      }
     1314      case TIME_FORMAT_MJD: {
     1315          double mjd = psTimeToMJD(time);
     1316          return psMetadataItemAllocF64(concept->name, concept->comment, mjd);
     1317      }
     1318      default: {
     1319          // If we've gotten this far, it's straight ISO.
     1320          psString dateTimeString = psTimeToISO(time); // String representation
     1321          psMetadataItem *item = psMetadataItemAllocStr(concept->name, concept->comment, dateTimeString);
     1322          psFree(dateTimeString);
     1323          return item;
     1324      }
     1325    }
     1326
     1327    return NULL;
    13361328}
    13371329

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/concepts/pmConceptsWrite.c

    r23594 r24244  
    6262        while ((cItem = psListGetAndIncrement(cIter))) {
    6363            if (cItem->type != PS_DATA_STRING) {
    64                 psWarning("psMetadataItem from list is of type %x instead of %x (PS_DATA_STRING) --- can't interpret.\n", cItem->type, PS_DATA_STRING);
     64                psWarning("psMetadataItem from list is of type %x instead of %x (PS_DATA_STRING) --- can't interpret.\n", cItem->type, PS_DATA_STRING);
    6565                psFree(cIter);
    6666                psFree(sIter);
    … …  
    606606        for (long i = 0; i < chips->n; i++) {
    607607            pmChip *chip = chips->data[i];  // Chip of interest
    608             if (chip && !chip->hdu) {
     608            if (chip) {
    609609                success &= pmConceptsWriteChip(chip, false, true, config);
    610610            }
    … …  
    622622    pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
    623623    bool success = conceptsWrite(&conceptsChip, fpa, chip, NULL, config, chip->concepts);
    624     if (propagateUp && !fpa->hdu) {
     624    if (propagateUp) {
    625625        psMetadata *conceptsFPA = pmConceptsSpecs(PM_FPA_LEVEL_FPA); // Concept specifications
    626626        success &= conceptsWrite(&conceptsFPA, fpa, chip, NULL, config, fpa->concepts);
    … …  
    630630        for (long i = 0; i < cells->n; i++) {
    631631            pmCell *cell = cells->data[i];  // Cell of interest
    632             if (cell && !cell->hdu) {
     632            if (cell) {
    633633                success &= pmConceptsWriteCell(cell, false, config);
    634634            }
    … …  
    649649    bool success = conceptsWrite(&conceptsCell, fpa, chip, cell, config, cell->concepts);
    650650    if (propagateUp) {
    651         if (!chip->hdu) {
    652             psMetadata *conceptsChip = pmConceptsSpecs(PM_FPA_LEVEL_CHIP); // Concept specifications
    653             success &= conceptsWrite(&conceptsChip, fpa, chip, cell, config, chip->concepts);
    654             if (!fpa->hdu) {
    655                 psMetadata *conceptsFPA = pmConceptsSpecs(PM_FPA_LEVEL_FPA); // Concept specifications
    656                 success &= conceptsWrite(&conceptsFPA, fpa, chip, cell, config, fpa->concepts);
    657             }
    658         }
     651        psMetadata *conceptsChip = pmConceptsSpecs(PM_FPA_LEVEL_CHIP); // Concept specifications
     652        success &= conceptsWrite(&conceptsChip, fpa, chip, cell, config, chip->concepts);
     653        psMetadata *conceptsFPA = pmConceptsSpecs(PM_FPA_LEVEL_FPA); // Concept specifications
     654        success &= conceptsWrite(&conceptsFPA, fpa, chip, cell, config, fpa->concepts);
    659655    }
    660656

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config

    • Property svn:ignore

      •  

        old new  
        77pmErrorCodes.c
        88pmErrorCodes.h
         9pmVersionDefinitions.h

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/Makefile.am

    r23351 r24244  
    11noinst_LTLIBRARIES = libpsmodulesconfig.la
    22
    3 # PSMODULES_VERSION=`if [ -e ../../VERSION ]; then cat ../../VERSION; else svnversion; fi`
    4 # PSMODULES_BRANCH=`if [ -e ../../BRANCH ]; then cat ../../BRANCH; else svn info | sed -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'; fi`
    5 # PSMODULES_SOURCE=`if [ -e ../../SOURCE ]; then cat ../../SOURCE; else svn info | sed -n -e 's/Repository UUID: // p'; fi`
     3if HAVE_SVNVERSION
     4PSMODULES_VERSION=`$(SVNVERSION) ../..`
     5else
     6PSMODULES_VERSION="UNKNOWN"
     7endif
     8
     9if HAVE_SVN
     10PSMODULES_BRANCH=`$(SVN) info ../.. | $(SED) -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'`
     11PSMODULES_SOURCE=`$(SVN) info | $(SED) -n -e 's/Repository UUID: // p'`
     12else
     13PSMODULES_BRANCH="UNKNOWN"
     14PSMODULES_SOURCE="UNKNOWN"
     15endif
    616
    717# Force recompilation of pmVersion.c, since it gets the version information
    8 # pmVersion.c: FORCE
    9 #       touch pmVersion.c
    10 # FORCE: ;
     18pmVersion.c: pmVersionDefinitions.h
     19pmVersionDefinitions.h: pmVersionDefinitions.h.in FORCE
     20        -$(RM) pmVersionDefinitions.h
     21        $(SED) -e "s|@PSMODULES_VERSION@|\"$(PSMODULES_VERSION)\"|" -e "s|@PSMODULES_BRANCH@|\"$(PSMODULES_BRANCH)\"|" -e "s|@PSMODULES_SOURCE@|\"$(PSMODULES_SOURCE)\"|" pmVersionDefinitions.h.in > pmVersionDefinitions.h
     22FORCE: ;
    1123
    12 libpsmodulesconfig_la_CPPFLAGS = $(SRCINC) $(PSMODULES_CFLAGS) -DPSMODULES_VERSION=$(SVN_VERSION) -DPSMODULES_BRANCH=$(SVN_BRANCH) -DPSMODULES_SOURCE=$(SVN_SOURCE)
     24libpsmodulesconfig_la_CPPFLAGS = $(SRCINC) $(PSMODULES_CFLAGS)
    1325libpsmodulesconfig_la_LDFLAGS  = -release $(PACKAGE_VERSION)
    1426libpsmodulesconfig_la_SOURCES  = \
    … …  
    3547
    3648# Error codes.
    37 BUILT_SOURCES = pmErrorCodes.h pmErrorCodes.c
    38 CLEANFILES = *~ pmErrorCodes.h pmErrorCodes.c
     49BUILT_SOURCES = pmErrorCodes.h pmErrorCodes.c pmVersionDefinitions.h
     50CLEANFILES = *~ pmErrorCodes.h pmErrorCodes.c pmVersionDefinitions.h
    3951
    4052pmErrorCodes.h : pmErrorCodes.dat pmErrorCodes.h.in

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfig.c

    r23594 r24244  
    551551                psWarning("Unable to resolve log destination: %s --- ignored", logDest);
    552552            } else {
     553                pmConfigRunFilenameAddWrite(config, "LOG", logDest);
    553554                config->logFD = psMessageDestination(resolved);
    554555            }
    … …  
    609610                psWarning("Unable to resolve trace destination: %s --- ignored", traceDest);
    610611            } else {
     612                pmConfigRunFilenameAddWrite(config, "TRACE", traceDest);
    611613                config->traceFD = psMessageDestination(resolved);
    612614            }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfigDump.c

    r23351 r24244  
    138138
    139139
    140 bool pmConfigDump(const pmConfig *config, const pmFPA *source, const char *filename)
     140bool pmConfigDump(const pmConfig *config, const char *filename)
    141141{
    142142    PS_ASSERT_PTR_NON_NULL(config, false);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfigDump.h

    r19399 r24244  
    3131///
    3232bool pmConfigDump(const pmConfig *config, ///< Configuration to dump
    33                   const pmFPA *source,    ///< Source FPA, defines the level for the file rule
    3433                  const char *filename    ///< Output file name
    3534    );

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfigRecipes.c

    r23351 r24244  
    8181            psTrace ("psModules.config", 3, "read recipes from camera config");
    8282        } else {
    83             psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config");
     83            psTrace ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config");
    8484        }
    8585    }
    … …  
    9797            psTrace ("psModules.config", 3, "read recipes from symbolic references");
    9898        } else {
    99             psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
     99            psTrace ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
    100100        }
    101101
    … …  
    105105            return false;
    106106        }
    107         psLogMsg ("psModules.config", PS_LOG_DETAIL, "merged camera recipes with system recipes");
     107        psTrace ("psModules.config", PS_LOG_DETAIL, "merged camera recipes with system recipes");
    108108
    109109        // load recipe-files specified on the command line
    … …  
    115115            psTrace ("psModules.config", 3, "read recipes from command-line arguments");
    116116        } else {
    117             psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
     117            psTrace ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
    118118        }
    119119
    … …  
    126126            psTrace ("psModules.config", 3, "read recipes from symbolic references");
    127127        } else {
    128             psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
     128            psTrace ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
    129129        }
    130130
    … …  
    137137            psTrace ("psModules.config", 3, "read recipes from command-line arguments");
    138138        } else {
    139             psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
     139            psTrace ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
    140140        }
    141141    }
    … …  
    348348            psFree(recipesIter);
    349349            return false;
    350         }
     350        }
    351351    }
    352352    psFree(recipesIter);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfigRun.c

    r23594 r24244  
    4141// Add a file to a nominated metadata in the RUN information
    4242static bool configRunFileAdd(pmConfig *config, // Configuration
    43                              const pmFPAfile *file, // File to add
     43                             const char *description, // Description of file
     44                             const char *name, // Name of file
    4445                             const char *target // Name of metadata to which to add
    4546                             )
    4647{
    4748    PS_ASSERT_PTR_NON_NULL(config, false);
    48     PS_ASSERT_PTR_NON_NULL(file, false);
     49    PS_ASSERT_STRING_NON_EMPTY(description, false);
     50    PS_ASSERT_STRING_NON_EMPTY(name, false);
     51    PS_ASSERT_STRING_NON_EMPTY(target, false);
    4952
    5053    psMetadata *run = configRun(config);// RUN information
    … …  
    5356    psAssert(files, "Require list of files");
    5457
    55     const char *name = file->name;      // Name of symbolic file
    56     const char *value = file->origname ? file->origname : file->filename; // The file (system) name
    57 
    5858    psString regex = NULL;              // Regular expression for iteration
    59     psStringAppend(&regex, "^%s$", name);
     59    psStringAppend(&regex, "^%s$", description);
    6060    psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex);
    6161    psFree(regex);
    … …  
    6363    while ((item = psMetadataGetAndIncrement(iter))) {
    6464        psAssert(item->type == PS_DATA_STRING, "We only put STRING types here.");
    65         if (strcmp(item->data.str, value) == 0) {
     65        if (strcmp(item->data.str, name) == 0) {
    6666            // It's already present
    6767            psFree(iter);
    … …  
    7171    psFree(iter);
    7272
    73     return psMetadataAddStr(files, PS_LIST_TAIL, name, PS_META_DUPLICATE_OK, NULL, value);
     73    return psMetadataAddStr(files, PS_LIST_TAIL, description, PS_META_DUPLICATE_OK, NULL, name);
    7474}
    7575
    … …  
    7979    PS_ASSERT_PTR_NON_NULL(file, false);
    8080
    81     return configRunFileAdd(config, file, "FILES.INPUT");
     81    return configRunFileAdd(config, file->name, file->origname ? file->origname : file->filename,
     82                            "FILES.INPUT");
     83}
     84
     85bool pmConfigRunFilenameAddRead(pmConfig *config, const char *description, const char *name)
     86{
     87    PS_ASSERT_PTR_NON_NULL(config, false);
     88    PS_ASSERT_STRING_NON_EMPTY(description, false);
     89    PS_ASSERT_STRING_NON_EMPTY(name, false);
     90
     91    return configRunFileAdd(config, description, name, "FILES.INPUT");
    8292}
    8393
    … …  
    8797    PS_ASSERT_PTR_NON_NULL(file, false);
    8898
    89     return configRunFileAdd(config, file, "FILES.OUTPUT");
     99    return configRunFileAdd(config, file->name, file->origname ? file->origname : file->filename,
     100                            "FILES.OUTPUT");
     101}
     102
     103bool pmConfigRunFilenameAddWrite(pmConfig *config, const char *description, const char *name)
     104{
     105    PS_ASSERT_PTR_NON_NULL(config, false);
     106    PS_ASSERT_STRING_NON_EMPTY(description, false);
     107    PS_ASSERT_STRING_NON_EMPTY(name, false);
     108
     109    return configRunFileAdd(config, description, name, "FILES.OUTPUT");
    90110}
    91111

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmConfigRun.h

    r23594 r24244  
    1212    );
    1313
     14/// Add a filename to the list of files read in the run-time information
     15bool pmConfigRunFilenameAddRead(
     16    pmConfig *config,                   ///< Configuration
     17    const char *description,            ///< Description of file
     18    const char *name                    ///< Name of file
     19    );
     20
    1421/// Add a file to the list of files written in the run-time information
    1522bool pmConfigRunFileAddWrite(
    1623    pmConfig *config,                   ///< Configuration
    1724    const pmFPAfile *file               ///< File to add
     25    );
     26
     27/// Add a filename to the list of files written in the run-time information
     28bool pmConfigRunFilenameAddWrite(
     29    pmConfig *config,                   ///< Configuration
     30    const char *description,            ///< Description of file
     31    const char *name                    ///< Name of file
    1832    );
    1933

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmErrorCodes.dat

    r13654 r24244  
    1111OBJECTS                 Problem in objects
    1212SKY                     Problem in sky
     13STAMPS                  Unable to select stamps for PSF-matching
    1314# these errors correspond to standard exit conditions
    1415ARGUMENTS               Incorrect arguments

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmErrorCodes.h.in

    r11300 r24244  
    99 */
    1010typedef enum {
    11     PM_ERR_BASE = 1200,
     11    PM_ERR_BASE = 1000,
    1212    PM_ERR_${ErrorCode},
    1313    PM_ERR_NERROR

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/config/pmVersion.c

    r23351 r24244  
    77#include <pslib.h>
    88#include "pmVersion.h"
     9#include "pmVersionDefinitions.h"
    910
    1011#ifndef PSMODULES_VERSION
    … …  
    1819#endif
    1920
    20 #define xstr(s) str(s)
    21 #define str(s) #s
    22 
    2321psString psModulesVersion(void)
    2422{
    2523    char *value = NULL;
    26     psStringAppend(&value, "%s@%s", xstr(PSMODULES_BRANCH), xstr(PSMODULES_VERSION));
     24    psStringAppend(&value, "%s@%s", PSMODULES_BRANCH, PSMODULES_VERSION);
    2725    return value;
    2826}
    … …  
    3028psString psModulesSource(void)
    3129{
    32     return psStringCopy(xstr(PSMODULES_SOURCE));
     30    return psStringCopy(PSMODULES_SOURCE);
    3331}
    3432

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmDetrendDB.c

    r23351 r24244  
    175175    if (!pipe) {
    176176        psError (PS_ERR_IO, false, "error calling command %s", line);
     177        psLogMsg ("psModules.detrend", PS_LOG_ERROR, "detselect command:\n %s\n", line);
    177178        goto failure;
    178179    }
    … …  
    188189        psError (PS_ERR_IO, false, "error running detselect");
    189190        psLogMsg ("psModules.detrend", PS_LOG_ERROR, "detselect command:\n %s\n", line);
     191        goto failure;
     192    }
     193
     194    if (!buffer->data || strlen(buffer->data) == 0) {
     195        psError(PS_ERR_IO, true, "Unable to find suitable detrend");
     196        psLogMsg("psModules.detrend", PS_LOG_ERROR, "detselect command:\n %s\n", line);
    190197        goto failure;
    191198    }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmFringeStats.c

    r23351 r24244  
    821821
    822822    // Solve the least-squares equation
    823     if (! psMatrixGJSolve(A, B)) {
     823    if (!psMatrixGJSolve(A, B)) {
    824824        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    825825        return false;
    … …  
    882882
    883883    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics
    884     psVectorStats(stats, diffs, NULL, mask, 1);
     884    if (!psVectorStats(stats, diffs, NULL, mask, 1)) {
     885        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     886        return 0;
     887    }
    885888    float middle = stats->sampleMedian; // The middle of the distribution
    886889    float thresh = rej * 0.74 * (stats->sampleUQ - stats->sampleLQ); // The rejection threshold
    … …  
    969972
    970973    // Get rid of the extreme outliers by assuming most of the points are somewhat clustered
    971     psVectorStats(median, science->f, NULL, NULL, 0);
     974    if (!psVectorStats(median, science->f, NULL, NULL, 0)) {
     975        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     976        return NULL;
     977    }
    972978    scale->coeff->data.F32[0] = median->sampleMedian;
    973979    for (int i = 0; i < fringes->n; i++) {
    974980        pmFringeStats *fringe = fringes->data[i]; // The fringe of interest
    975         psVectorStats(median, fringe->f, NULL, NULL, 0);
     981        if (!psVectorStats(median, fringe->f, NULL, NULL, 0)) {
     982            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     983            return NULL;
     984        }
    976985        scale->coeff->data.F32[0] -= median->sampleMedian;
    977986        scale->coeff->data.F32[i] = 0.0;

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmMaskBadPixels.c

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmMaskBadPixels.h

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmOverscan.c

    r21183 r24244  
    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) {
    … …  
    255258    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
    256259                     comment, "");
     260    psTrace ("psModules.detrend", 4, "%s\n", comment);
    257261    psFree(comment);
    258262
    … …  
    274278        psFree(iter);
    275279
    276         (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        }
    277284        psFree(pixels);
    278285        double reduced = psStatsGetValue(stats, statistic); // Result of statistics
    … …  
    344351        }
    345352
     353        // generate stats of overscan vector for header
     354        {
     355          psString comment = NULL;    // Comment to add
     356          psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     357          if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
     358              psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     359              return false;
     360          }
     361          psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
     362          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     363          psFree(comment);
     364
     365          // write metadata header value
     366          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
     367          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
     368          psFree (vectorStats);
     369        }
     370
    346371        // Subtract row by row
    347372        for (int i = 0; i < image->numRows; i++) {
    … …  
    392417        }
    393418
    394         // Subtract column by column
    395         for (int i = 0; i < image->numCols; i++) {
    396             for (int j = 0; j < image->numRows; j++) {
     419        // generate stats of overscan vector for header
     420        {
     421          psString comment = NULL;    // Comment to add
     422          psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     423          if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
     424              psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     425              return false;
     426          }
     427          psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
     428          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     429          psFree(comment);
     430
     431          // write metadata header value
     432          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
     433          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
     434          psFree (vectorStats);
     435        }
     436
     437        // Subtract column by column
     438        for (int j = 0; j < image->numRows; j++) {
     439          for (int i = 0; i < image->numCols; i++) {
    397440                image->data.F32[j][i] -= reduced->data.F32[i];
    398441            }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmRemnance.c

    r23351 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmShutterCorrection.c

    r23594 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/detrend/pmSkySubtract.c

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/extras/pmVisual.c

    r23351 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmImageCombine.c

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmPSFEnvelope.c

    r23351 r24244  
    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);
    … …  
    265299
    266300        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         }
     301            // Can't do anything about it; limp along as best we can
     302            psErrorClear();
     303            continue;
     304        }
     305        numMoments++;
     306    }
     307
     308    if (numMoments == 0) {
     309        psError(PS_ERR_UNKNOWN, true, "Unable to measure moments for sources.");
     310        psFree(fakes);
     311        psFree(readout);
     312        return NULL;
    272313    }
    273314
    … …  
    286327    options->psfFieldXo = 0;
    287328    options->psfFieldYo = 0;
     329    options->chiFluxTrend = false;      // All sources have similar flux, so fitting a trend often fails
    288330
    289331    pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmReadoutCombine.c

    r21363 r24244  
    384384            // Combination
    385385            if (!psVectorStats(stats, pixels, errors, mask, 1)) {
    386                 // Can't do much about it, but it's not worth worrying about
    387                 psErrorClear();
     386                psError(PS_ERR_UNKNOWN, false, "error in pixel stats");
     387                return false;
     388            }
     389           
     390            outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
     391
     392            if (isnan(outputImage[yOut][xOut])) {
    388393                outputImage[yOut][xOut] = NAN;
    389394                outputMask[yOut][xOut] = params->blank;
     395                sigma->data.F32[yOut][xOut] = NAN;
    390396                if (params->variances) {
    391397                    outputVariance[yOut][xOut] = NAN;
    392398                }
    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             }
     399                continue;
     400            }
     401            outputMask[yOut][xOut] = 0;
     402            sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev);
     403            if (params->variances) {
     404                float stdev = psStatsGetValue(stats, combineStdev);
     405                outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance
     406                // XXXX this is not the correct formal error.
     407                // also, the weighted mean is not obviously the correct thing here
     408            }
    405409        }
    406410    }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmStack.c

    r23594 r24244  
    3434
    3535//#define TESTING                         // Enable test output
    36 //#define TEST_X 2318                     // x coordinate to examine
    37 //#define TEST_Y 2306                     // y coordinate to examine
     36//#define TEST_X 1085                     // x coordinate to examine
     37//#define TEST_Y 3371                     // y coordinate to examine
    3838
    3939
    … …  
    265265
    266266    int numBad = frac * numGood + 0.5;  // Number of bad values
    267     int low = numBad / 2, high = low + numGood; // Indices (modulo masked pixels) delimiting range of interest
     267    int low = numBad / 2, high = low + numGood - numBad; // Indices (modulo masked pixels)
    268268
    269269    sortBuffer = psVectorSortIndex(sortBuffer, values);
    … …  
    536536#ifdef TESTING
    537537                  if (x == TEST_X && y == TEST_Y) {
    538                       fprintf(stderr, "Rejection limit: %f\n", limit);
     538                      fprintf(stderr, "Rejecting without variance; rejection limit: %f\n", limit);
    539539                  }
    540540#endif
    541541              } else {
     542#ifdef TESTING
     543                  if (x == TEST_X && y == TEST_Y) {
     544                      fprintf(stderr, "Rejecting with variance...\n");
     545                  }
     546#endif
    542547                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
    543548              }

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtraction.c

    r21470 r24244  
    822822    psFree(mask);
    823823
     824    // XXX raise an error?
     825    if (isnan(stats->sampleMean)) {
     826        psFree(stats);
     827        return -1;
     828    }
     829
    824830    double mean, rms;                 // Mean and RMS of deviations
    825831    if (numStamps < MIN_SAMPLE_STATS) {
    … …  
    849855    float limit = sigmaRej * rms; // Limit on maximum deviation
    850856    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
     857
     858
     859    psString ds9name = NULL;            // Filename for ds9 region file
     860    static int ds9num = 0;              // File number for ds9 region file
     861    psStringAppend(&ds9name, "stamps_reject_%d.ds9", ds9num);
     862    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "rejected stamps");
     863    psFree(ds9name);
     864    ds9num++;
    851865
    852866    int numRejected = 0;                // Number of stamps rejected
    … …  
    872886                    }
    873887                }
     888                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    874889
    875890                // Set stamp for replacement
    … …  
    897912                numGood++;
    898913                newMean += deviations->data.F32[i];
     914                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    899915            }
    900916        }
    901917    }
    902918    newMean /= numGood;
     919
     920    if (ds9) {
     921        fclose(ds9);
     922    }
    903923
    904924    if (numRejected > 0) {
    … …  
    13851405    }
    13861406
     1407    // Data exists on the outputs now
     1408    if (out1) {
     1409        out1->data_exists = true;
     1410        if (out1->parent) {
     1411            out1->parent->data_exists = out1->parent->parent->data_exists = true;
     1412        }
     1413    }
     1414    if (out2) {
     1415        out2->data_exists = true;
     1416        if (out2->parent) {
     1417            out2->parent->data_exists = out2->parent->parent->data_exists = true;
     1418        }
     1419    }
     1420
     1421
    13871422    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve image: %f sec",
    13881423             psTimerClear("pmSubtractionConvolve"));

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionAnalysis.c

    r21422 r24244  
    201201        psFree(image);
    202202
    203         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX,
    204                          PS_META_DUPLICATE_OK, "Maximum deconvolution fraction", max);
     203        psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
     204        if (item) {
     205            item->data.F32 = PS_MAX(item->data.F32, max);
     206        } else {
     207            psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
     208                             "Maximum deconvolution fraction", max);
     209        }
    205210    }
    206211

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionEquation.c

    r21422 r24244  
    749749    }
    750750
     751    psString ds9name = NULL;            // Filename for ds9 region file
     752    static int ds9num = 0;              // File number for ds9 region file
     753    psStringAppend(&ds9name, "stamps_solution_%d.ds9", ds9num);
     754    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "solution stamps");
     755    psFree(ds9name);
     756    ds9num++;
     757
    751758    if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
    752759        // Accumulate the least-squares matricies and vectors
    … …  
    761768                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
    762769                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
     770                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    763771                numStamps++;
     772            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
     773                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    764774            }
    765775        }
    … …  
    788798
    789799        psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    790         psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
     800        psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    791801        psFree(sumMatrix);
    792802        if (!luMatrix) {
    … …  
    797807            return NULL;
    798808        }
    799         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
     809        kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
    800810
    801811        psFree(sumVector);
    … …  
    830840                (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
    831841                (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
     842                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    832843                numStamps++;
    833844            }
    … …  
    9901001        // XXXXX Free temporary matrices and vectors
    9911002
     1003    }
     1004
     1005    if (ds9) {
     1006        fclose(ds9);
    9921007    }
    9931008

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionMask.c

    r21183 r24244  
    3838
    3939psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psImageMaskType maskVal,
    40                            int size, int footprint, float badFrac, bool useFFT)
     40                           int size, int footprint, float badFrac, pmSubtractionMode mode)
    4141{
    4242    PS_ASSERT_IMAGE_NON_NULL(mask1, NULL);
    … …  
    114114    }
    115115
    116     // XXX Could do something smarter here --- we will get images that are predominantly masked (where the
    117     // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad
    118     // pixel is wasting time.  As a first cut, I've put in a check on the fraction of bad pixels, but we could
    119     // imagine looking for the edge of big regions and convolving just at the edge.  As a second cut, allow
    120     // use of FFT convolution.
    121 
    122116    for (int y = 0; y < numRows; y++) {
    123117        for (int x = 0; x < numCols; x++) {
    … …  
    131125    }
    132126
    133     // Block out the entire stamp footprint around bad input pixels.
    134 
    135127    // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad
    136     // reference pixel (within 'size').  Then we want to block out with the FOOTPRINT mask everything within a
     128    // reference pixel (within 'size').  Then we want to block out with the REJ mask everything within a
    137129    // footprint's distance of those (within 'footprint').
    138130
    139131    bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolve
    140132
    141     if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1,
    142                             PM_SUBTRACTION_MASK_CONVOLVE_1,
     133    // Pixels that will be bad (or poor) if we convolve with a bad reference pixel
     134    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, PM_SUBTRACTION_MASK_CONVOLVE_1,
    143135                             -size, size, -size, size)) {
    144136        psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1.");
    … …  
    146138        return NULL;
    147139    }
    148     if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2,
    149                              PM_SUBTRACTION_MASK_CONVOLVE_2,
     140    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, PM_SUBTRACTION_MASK_CONVOLVE_2,
    150141                             -size, size, -size, size)) {
    151142        psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2.");
     143        psFree(mask);
     144        return NULL;
     145    }
     146
     147    // Pixels that should not be chosen as stamps
     148    psImageMaskType maskRej = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2 |
     149        PM_SUBTRACTION_MASK_BORDER;     // Mask value for rejection
     150    switch (mode) {
     151      case PM_SUBTRACTION_MODE_1:
     152        maskRej |= PM_SUBTRACTION_MASK_CONVOLVE_1;
     153        break;
     154      case PM_SUBTRACTION_MODE_2:
     155        maskRej |= PM_SUBTRACTION_MASK_CONVOLVE_2;
     156        break;
     157      case PM_SUBTRACTION_MODE_UNSURE:
     158      case PM_SUBTRACTION_MODE_DUAL:
     159        maskRej |= PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
     160        break;
     161      default:
     162        psAbort("Unsupported subtraction mode: %x", mode);
     163    }
     164    if (!psImageConvolveMask(mask, mask, maskRej, PM_SUBTRACTION_MASK_REJ,
     165                             -footprint, footprint, -footprint, footprint)) {
     166        psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels.");
    152167        psFree(mask);
    153168        return NULL;

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionMask.h

    r21183 r24244  
    1111                           int footprint, ///< Half-size of the kernel footprint
    1212                           float badFrac, ///< Maximum fraction of bad input pixels to accept
    13                            bool useFFT  ///< Use FFT to do convolution?
     13                           pmSubtractionMode mode  ///< Subtraction mode
    1414    );
    1515

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionMatch.c

    r23351 r24244  
    1818#include "pmSubtractionMask.h"
    1919#include "pmSubtractionThreads.h"
     20#include "pmSubtractionVisual.h"
     21#include "pmErrorCodes.h"
     22
    2023#include "pmSubtractionMatch.h"
    21 #include "pmSubtractionVisual.h"
    2224
    2325#define BG_STAT PS_STAT_ROBUST_MEDIAN   // Statistic to use for background
    … …  
    6264                      psImage *variance,  // Variance map
    6365                      const psRegion *region, // Region of interest, or NULL
    64                       float threshold,  // Threshold for stamp finding
     66                      float thresh1,  // Threshold for stamp finding on readout 1
     67                      float thresh2,  // Threshold for stamp finding on readout 2
    6568                      float stampSpacing, // Spacing between stamps
    6669                      int size,         // Kernel half-size
    … …  
    7073{
    7174    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    72     *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, footprint,
    73                                       stampSpacing, mode);
     75
     76    psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest
     77
     78    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
     79                                      size, footprint, stampSpacing, mode);
    7480    if (!*stamps) {
    75         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
     81        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    7682        return false;
    7783    }
    … …  
    238244
    239245    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
    240                                          badFrac, useFFT); // Subtraction mask
     246                                         badFrac, mode); // Subtraction mask
    241247    if (!subMask) {
    242248        psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
    … …  
    371377
    372378    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    373                                 badFrac, useFFT);
     379                                badFrac, subMode);
    374380    if (!subMask) {
    375381        psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
    … …  
    390396    }
    391397
     398    float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
    392399    {
    393400        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
    394         if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
    395             psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
    396             psFree(bg);
    397             psFree(rng);
    398             goto MATCH_ERROR;
    399         }
    400         threshold = bg->robustMedian + threshold * bg->robustStdev;
     401        if (ro1) {
     402            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
     403                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
     404                psFree(bg);
     405                goto MATCH_ERROR;
     406            }
     407            stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
     408        }
     409        if (ro2) {
     410            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
     411                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
     412                psFree(bg);
     413                goto MATCH_ERROR;
     414            }
     415            stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
     416        }
    401417        psFree(bg);
    402418    }
    … …  
    416432            }
    417433
    418             if (sources) {
    419                 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, footprint,
    420                                                           stampSpacing, subMode);
    421             } else if (stampsName && strlen(stampsName) > 0) {
    422                 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint,
    423                                                         stampSpacing, subMode);
     434            if (stampsName && strlen(stampsName) > 0) {
     435                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
     436                                                        footprint, stampSpacing, subMode);
     437            } else if (sources) {
     438                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
     439                                                           footprint, stampSpacing, subMode);
    424440            }
    425441
    426442            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
    427443            // doesn't matter.
    428             if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, threshold, stampSpacing,
    429                            size, footprint, subMode)) {
     444            if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
     445                           stampSpacing, size, footprint, subMode)) {
    430446                goto MATCH_ERROR;
    431447            }
    … …  
    488504                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    489505
    490                 if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, threshold, stampSpacing,
    491                                size, footprint, subMode)) {
     506                if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
     507                               stampSpacing, size, footprint, subMode)) {
    492508                    goto MATCH_ERROR;
    493509                }
    … …  
    565581            psFree(kernels);
    566582            kernels = NULL;
    567 
    568             // There is data in the readout now
    569             if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    570                 conv1->data_exists = true;
    571                 if (conv1->parent) {
    572                     conv1->parent->data_exists = true;
    573                     conv1->parent->parent->data_exists = true;
    574                 }
    575             }
    576             if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    577                 conv2->data_exists = true;
    578                 if (conv2->parent) {
    579                     conv2->parent->data_exists = true;
    580                     conv2->parent->parent->data_exists = true;
    581                 }
    582             }
    583583        }
    584584    }
    … …  
    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/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionStamps.c

    r21363 r24244  
    2020#include "pmModel.h"
    2121#include "pmSource.h"
    22 
     22#include "pmErrorCodes.h"
    2323
    2424#include "pmSubtraction.h"
    … …  
    7878}
    7979
    80 // Is this position unmasked?
    81 static bool checkStampMask(int x, int y, // Coordinates of stamp
    82                            const psImage *mask, // Mask
    83                            pmSubtractionMode mode, // Mode for subtraction
    84                            int footprint // Footprint to check for Bad Stuff
    85                            )
    86 {
    87     if (!mask) {
    88         return true;
    89     }
    90 
    91     bool clean = true;                  // Is the footprint clean?
    92     int numCols = mask->numCols, numRows = mask->numRows; // Size of image
    93 
    94     // Check the footprint bounds
    95     if (x < footprint || x >= numCols - footprint || y < footprint || y >= numRows - footprint) {
    96         clean = false;
    97     }
    98 
    99     // Determine mask value
    100     psImageMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
    101     switch (mode) {
    102       case PM_SUBTRACTION_MODE_1:
    103         maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1;
    104         break;
    105       case PM_SUBTRACTION_MODE_2:
    106         maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_2;
    107         break;
    108       case PM_SUBTRACTION_MODE_UNSURE:
    109       case PM_SUBTRACTION_MODE_DUAL:
    110         maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
    111         break;
    112       default:
    113         psAbort("Unsupported subtraction mode: %x", mode);
    114     }
    115 
    116     // Check the immediate pixel
    117     if (clean && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (maskVal | PM_SUBTRACTION_MASK_REJ))) {
    118         clean = false;
    119     }
    120 
    121     int xMin = PS_MAX(x - footprint, 0), xMax = PS_MIN(x + footprint, numCols - 1); // Bounds in x
    122     int yMin = PS_MAX(y - footprint, 0), yMax = PS_MIN(y + footprint, numRows - 1); // Bounds in y
    123 
    124     // Check the footprint
    125     if (clean) {
    126         for (int j = yMin; j <= yMax; j++) {
    127             for (int i = xMin; i <= xMax; i++) {
    128                 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal) {
    129                     clean = false;
    130                     goto CHECK_STAMP_MASK_DONE;
    131                 }
     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;
    132106            }
    133         }
    134     }
    135 CHECK_STAMP_MASK_DONE:
    136 
    137     if (!clean) {
    138         // Mask the footprint, so we don't select something near it again
    139         for (int j = yMin; j <= yMax; j++) {
    140             for (int i = xMin; i <= xMax; i++) {
    141                 mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= PM_SUBTRACTION_MASK_REJ;
     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;
    142111            }
    143         }
    144     }
    145 
    146     return clean;
    147 }
     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;
     126}
     127
     128//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     129// Functions for generating ds9 region files
     130//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     131
     132static bool ds9regions = false;         // Save ds9 region files?
     133
     134void pmSubtractionRegions(bool state)
     135{
     136    ds9regions = state;
     137}
     138
     139FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, const char *filename,
     140                              const char *description)
     141{
     142    if (!ds9regions || !stamps || !filename) {
     143        return NULL;
     144    }
     145
     146    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Writing %s to ds9 region file: %s",
     147             description, filename);
     148
     149    FILE *file = fopen(filename, "w");
     150
     151    // Outline the stamps
     152    for (int i = 0; i < stamps->num; i++) {
     153        psRegion *region = stamps->regions->data[i]; // Region of interest
     154        float xCentre = 0.5 * (region->x0 + region->x1), yCentre = 0.5 * (region->y0 + region->y1);
     155        fprintf(file, "image;box(%f,%f,%f,%f,0) # color=blue\nimage;text(%f,%f,{%d}) # color=blue\n",
     156                xCentre, yCentre, region->x1 - region->x0, region->y1 - region->y0,
     157                xCentre, yCentre, i);
     158    }
     159
     160    return file;
     161}
     162
     163void pmSubtractionStampPrint(FILE *ds9, float x, float y, float size, const char *color)
     164{
     165    if (!ds9regions || !ds9) {
     166        return;
     167    }
     168    fprintf(ds9, "image;circle(%f,%f,%f)", x, y, size);
     169    if (color && strlen(color) > 0) {
     170        fprintf(ds9, " # color=%s", color);
     171    }
     172    fprintf(ds9, "\n");
     173    return;
     174}
     175
    148176
    149177//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    225253
    226254
    227 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
    228                                                 const psImage *subMask, const psRegion *region,
    229                                                 float threshold, int footprint, float spacing,
     255pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1,
     256                                                const psImage *image2, const psImage *subMask,
     257                                                const psRegion *region, float thresh1, float thresh2,
     258                                                int size, int footprint, float spacing,
    230259                                                pmSubtractionMode mode)
    231260{
    232     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    233     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
     261    if (!image1 && !image2) {
     262        psError(PS_ERR_UNEXPECTED_NULL, true, "Must specify an image");
     263        return NULL;
     264    }
     265    int numCols = 0, numRows = 0;       // Size of images
     266    if (image1) {
     267        PS_ASSERT_IMAGE_NON_NULL(image1, NULL);
     268        PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, NULL);
     269        if (subMask) {
     270            PS_ASSERT_IMAGES_SIZE_EQUAL(image1, subMask, NULL);
     271        }
     272        numCols = image1->numCols;
     273        numRows = image1->numRows;
     274    }
     275    if (image2) {
     276        PS_ASSERT_IMAGE_NON_NULL(image2, NULL);
     277        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, NULL);
     278        if (subMask) {
     279            PS_ASSERT_IMAGES_SIZE_EQUAL(image2, subMask, NULL);
     280        }
     281        numCols = image2->numCols;
     282        numRows = image2->numRows;
     283    }
     284    if (image1 && image2) {
     285        PS_ASSERT_IMAGES_SIZE_EQUAL(image1, image2, NULL);
     286    }
    234287    if (subMask) {
    235288        PS_ASSERT_IMAGE_NON_NULL(subMask, NULL);
    236         PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);
    237289        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, NULL);
    238     }
     290        PS_ASSERT_IMAGE_SIZE(subMask, numCols, numRows, NULL);
     291    }
     292    PS_ASSERT_INT_POSITIVE(size, NULL);
    239293    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
    240294    PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
    … …  
    246300            return false;
    247301        }
    248         if (region->x0 < 0 || region->x1 > image->numCols ||
    249             region->y0 < 0 || region->y1 > image->numRows) {
     302        if (region->x0 < 0 || region->x1 > numCols ||
     303            region->y0 < 0 || region->y1 > numRows) {
    250304            psString string = psRegionToString(*region);
    251305            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
    252                     string, image->numCols, image->numRows);
     306                    string, numCols, numRows);
    253307            psFree(string);
    254308            return false;
    … …  
    256310    }
    257311
    258     int numRows = image->numRows, numCols = image->numCols; // Size of image
     312    int border = size + footprint;      // Border size
    259313
    260314    if (!stamps) {
    … …  
    281335            numSearch++;
    282336
    283             float xStamp = 0, yStamp = 0;  // Coordinates of stamp
    284             float fluxStamp = NAN;          // Flux of stamp
    285             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?
    286340
    287341            // A couple different ways of finding stamps:
    … …  
    291345                psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
    292346
    293                 // Take stamp off the top of the (sorted) list
    294                 if (xList->n > 0) {
    295                     int index = xList->n - 1; // Index of new stamp
    296                     xStamp = xList->data.F32[index];
    297                     yStamp = yList->data.F32[index];
    298                     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
    299350
    300351                    // Chop off the top of the list
    301                     xList->n = index;
    302                     yList->n = index;
    303                     fluxList->n = index;
    304 
    305                     goodStamp = (fluxStamp > threshold) ? true : false;
    306                 } else {
    307                     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);
    308369                }
    309370            } else {
    310371                // Use a simple method of automatically finding stamps --- take the highest pixel in the
    311372                // subregion
    312                 fluxStamp = threshold;
    313373                psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
    314                 for (int y = subRegion->y0; y <= subRegion->y1; y++) {
    315                     for (int x = subRegion->x0; x <= subRegion->x1; x++) {
    316                         if (checkStampMask(x, y, subMask, mode, footprint) &&
    317                             image->data.F32[y][x] > fluxStamp) {
    318                             fluxStamp = image->data.F32[y][x];
    319                             xStamp = x;
    320                             yStamp = y;
    321                             goodStamp = true;
    322                         }
    323                     }
    324                 }
     374
     375                goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
     376                                        subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1,
     377                                        numCols, numRows, border);
    325378            }
    326379
    … …  
    355408    if (numGood == 0 && numFound == 0) {
    356409        // No good stamps
     410        psError(PM_ERR_STAMPS, true, "Unable to find suitable stamps");
    357411        psFree(stamps);
    358412        return NULL;
    … …  
    363417
    364418
     419
    365420pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y,
    366421                                               const psImage *image, const psImage *subMask,
    367                                                const psRegion *region, int footprint, float spacing,
    368                                                pmSubtractionMode mode)
     422                                               const psRegion *region, int size, int footprint,
     423                                               float spacing, pmSubtractionMode mode)
    369424
    370425{
    … …  
    378433        PS_ASSERT_IMAGE_NON_NULL(subMask, NULL);
    379434        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, NULL);
    380         if (image) {
    381             PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    382             PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);
    383         }
    384     }
     435        PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);
     436    }
     437    PS_ASSERT_INT_POSITIVE(size, NULL);
     438    PS_ASSERT_INT_POSITIVE(footprint, NULL);
    385439    PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
    386440
    … …  
    389443                                                                 region, footprint, spacing); // Stamp list
    390444    int numStamps = stamps->num;        // Number of stamps
     445
     446    psString ds9name = NULL;            // Filename for ds9 region file
     447    static int ds9num = 0;              // File number for ds9 region file
     448    psStringAppend(&ds9name, "stamps_all_%d.ds9", ds9num);
     449    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "all stamps"); // ds9 region file
     450    ds9num++;
    391451
    392452    // Initialise the lists
    … …  
    408468            psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region",
    409469                    xPix, yPix);
    410             continue;
    411         }
    412         if (!checkStampMask(xPix, yPix, subMask, mode, footprint)) {
    413             // Not a good stamp
    414             psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",
    415                     xPix, yPix);
     470            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red");
    416471            continue;
    417472        }
    … …  
    437492                psTrace("psModules.imcombine", 9, "Putting input stamp (%d,%d) into subregion %d",
    438493                        xPix, yPix, j);
     494                pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "green");
    439495            }
    440496        }
    … …  
    443499            psTrace("psModules.imcombine", 9, "Unable to find subregion for stamp (%d,%d)",
    444500                    xPix, yPix);
    445         }
     501            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "yellow");
     502        }
     503    }
     504
     505    if (ds9) {
     506        fclose(ds9);
    446507    }
    447508
    … …  
    474535    }
    475536
    476 
    477537    return stamps;
    478538}
    … …  
    544604}
    545605
    546 #if 0
    547 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize)
    548 {
    549     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    550     PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false);
    551     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    552 
    553     int size = kernelSize + stamps->footprint; // Size of postage stamps
    554     int num = stamps->num;              // Number of stamps
    555     float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma
    556 
    557     for (int i = 0; i < num; i++) {
    558         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    559         if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) {
    560             continue;
    561         }
    562 
    563         float x = stamp->x, y = stamp->y; // Coordinates of stamp
    564         float flux = stamp->flux; // Flux of star
    565         if (!isfinite(flux)) {
    566             psWarning("Unable to generate PSF for stamp %d --- bad flux.", i);
    567             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    568             continue;
    569         }
    570 
    571         float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame
    572         float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
    573 
    574         psFree(stamp->image2);
    575         stamp->image2 = psKernelAlloc(-size, size, -size, size);
    576         psKernel *target = stamp->image2; // Target stamp
    577 
    578         // Put in a Waussian, just for fun!
    579         for (int v = -size; v <= size; v++) {
    580             for (int u = -size; u <= size; u++) {
    581                 float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
    582                 target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
    583             }
    584         }
    585 
    586     }
    587 
    588     return true;
    589 }
    590 #endif
    591 
    592606pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
    593607                                                          const psImage *subMask, const psRegion *region,
    594                                                           int footprint, float spacing,
     608                                                          int size, int footprint, float spacing,
    595609                                                          pmSubtractionMode mode)
    596610{
    … …  
    621635    y->n = numOut;
    622636
    623     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region,
     637    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
    624638                                                            footprint, spacing, mode); // Stamps to return
    625639    psFree(x);
    … …  
    636650pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
    637651                                                       const psImage *subMask, const psRegion *region,
    638                                                        int footprint, float spacing, pmSubtractionMode mode)
     652                                                       int size, int footprint, float spacing,
     653                                                       pmSubtractionMode mode)
    639654{
    640655    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
    … …  
    652667    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    653668
    654     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, footprint,
     669    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
    655670                                                            spacing, mode);
    656671    psFree(data);

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionStamps.h

    r21363 r24244  
    7474/// Find stamps on an image
    7575pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
    76                                                 const psImage *image, ///< Image for which to find stamps
     76                                                const psImage *image1, ///< Image for which to find stamps
     77                                                const psImage *image2, ///< Image for which to find stamps
    7778                                                const psImage *mask, ///< Mask, or NULL
    7879                                                const psRegion *region, ///< Region to search, or NULL
    79                                                 float threshold, ///< Threshold for stamps in the image
     80                                                float thresh1, ///< Threshold for stamps in image 1
     81                                                float thresh2, ///< Threshold for stamps in image 2
     82                                                int size, ///< Kernel half-size
    8083                                                int footprint, ///< Half-size for stamps
    8184                                                float spacing, ///< Rough spacing for stamps
    … …  
    8992                                               const psImage *mask, ///< Mask, or NULL
    9093                                               const psRegion *region, ///< Region to search, or NULL
     94                                               int size, ///< Kernel half-size
    9195                                               int footprint, ///< Half-size for stamps
    9296                                               float spacing, ///< Rough spacing for stamps
    … …  
    100104    const psImage *subMask,             ///< Mask, or NULL
    101105    const psRegion *region,             ///< Region to search, or NULL
     106    int size,                           ///< Kernel half-size
    102107    int footprint,                      ///< Half-size for stamps
    103108    float spacing,                      ///< Rough spacing for stamps
    … …  
    111116    const psImage *subMask,             ///< Mask, or NULL
    112117    const psRegion *region,             ///< Region to search, or NULL
     118    int size,                           ///< Kernel half-size
    113119    int footprint,                      ///< Half-size for stamps
    114120    float spacing,                      ///< Rough spacing for stamps
    … …  
    124130    );
    125131
     132
     133/// Turn on/off generation of ds9 region files
     134///
     135/// Intended for debugging
     136void pmSubtractionRegions(bool state    ///< Generate ds9 region files?
     137    );
     138
     139/// Open a file for ds9 regions
     140///
     141/// Intended for debugging
     142FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, ///< List of stamps, for outlines
     143                              const char *filename, ///< Filename to write
     144                              const char *description ///< Description of file
     145    );
     146
     147/// Print a stamp position to ds9 region file
     148///
     149/// Intended for debugging
     150void pmSubtractionStampPrint(FILE *ds9, ///< ds9 region file
     151                             float x, float y, ///< Position of stamp
     152                             float size,///< Size of circle
     153                             const char *color ///< Colour
     154    );
     155
    126156#endif

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/Makefile.am

    r23351 r24244  
    3838        pmSourceIO_PS1_CAL_0.c \
    3939        pmSourceIO_CMF_PS1_V1.c \
     40        pmSourceIO_MatchedRefs.c \
    4041        pmSourcePlots.c \
    4142        pmSourcePlotPSFModel.c \

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/models/pmModel_PS1_V1.c

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

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/models/pmModel_SGAUSS.c

    r15834 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmGrowthCurveGenerate.c

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmPSF.c

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

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmPSF.h

    r23594 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmPSFtry.c

    r21183 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSource.c

    r23351 r24244  
    479479        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    480480
    481         psVectorStats (stats, tmpSx, NULL, NULL, 0);
     481        if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) {
     482            psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
     483            return (emptyClump);
     484        }
    482485        psfClump.X  = stats->clippedMean;
    483486        psfClump.dX = stats->clippedStdev;
    484487
    485         psVectorStats (stats, tmpSy, NULL, NULL, 0);
     488        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
     489            psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
     490            return (emptyClump);
     491        }
    486492        psfClump.Y  = stats->clippedMean;
    487493        psfClump.dY = stats->clippedStdev;
    … …  
    636642
    637643        if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) {
    638             // Don't care about this error
    639             psErrorClear();
     644            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     645            psFree (stats);
     646            psFree (starsn_peaks);
     647            return false;
    640648        }
    641649        psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats->min, stats->max);
    … …  
    648656        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    649657        if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) {
    650             // Don't care about this error
    651             psErrorClear();
     658            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     659            psFree (stats);
     660            psFree (starsn_peaks);
     661            return false;
    652662        }
    653663        psLogMsg ("psModules.objects", 3, "SN range (peaks)  : %f - %f (%ld)\n",

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceFitModel.c

    r23351 r24244  
    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/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceIO.c

    r22699 r24244  
    562562
    563563    // not needed if only one chip
    564     if (file->fpa->chips->n == 1) return true;
    565 
     564    if (file->fpa->chips->n == 1) {
     565        pmSourceWriteMatchedRefs (file->fits, file->fpa, config);
     566        return true;
     567    }
    566568
    567569    // find the FPA phu
    … …  
    667669    psFree (outhead);
    668670
     671    pmSourceWriteMatchedRefs (file->fits, file->fpa, config);
    669672    return true;
    670673}
    … …  
    10531056}
    10541057
     1058   

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceIO.h

    r21516 r24244  
    7575bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos);
    7676
     77bool pmSourceWriteMatchedRefs (psFits *fits, pmFPA *fpa, pmConfig *config);
     78
    7779/// @}
    7880# endif /* PM_SOURCE_IO_H */

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/objects/pmSourceMatch.c

    r23351 r24244  
    473473        return -1;
    474474    }
     475    // XXX handle this case better:
     476    if (isnan(stats->clippedMean))  {
     477        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
     478        psFree(stats);
     479        return -1;
     480    }
    475481
    476482    float thresh = stats->clippedMean + photoLevel * stats->clippedStdev; // Threshold for clouds

  • branches/cnb_branches/cnb_branch_20090301/psModules/src/psmodules.h

    r23351 r24244  
    2828#include <pmConfigMask.h>
    2929#include <pmConfigDump.h>
     30#include <pmConfigRun.h>
    3031#include <pmVersion.h>
    3132
Note: See TracChangeset for help on using the changeset viewer.