IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 24, 2010, 2:59:09 PM (16 years ago)
Author:
Paul Price
Message:

Merging trunk in advance of reintegrating into trunk.

Location:
branches/pap
Files:
1 deleted
21 edited

Legend:

Unmodified
Added
Removed

  • branches/pap

  • branches/pap/psModules

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

    r27862 r28484  
    378378    // XXX make it optional to write out CDi_j terms, or other versions
    379379    // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
    380     if (!wcs->wcsCDkeys) { 
     380    if (!wcs->wcsCDkeys) {
    381381
    382382      double cdelt1 = wcs->cdelt1;
    … …  
    384384      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
    385385      psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
    386      
     386
    387387      // test the PC00i00j varient:
    388388      psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
    … …  
    390390      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
    391391      psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
    392      
     392
    393393      // Elixir-style polynomial terms
    394394      // XXX currently, Elixir/DVO cannot accept mixed orders
    … …  
    398398      if (fitOrder > 1) {
    399399        for (int i = 0; i <= fitOrder; i++) {
    400           for (int j = 0; j <= fitOrder; j++) {
    401             if (i + j < 2)
    402               continue;
    403             if (i + j > fitOrder)
    404               continue;
    405             sprintf (name, "PCA1X%1dY%1d", i, j);
    406             psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    407             sprintf (name, "PCA2X%1dY%1d", i, j);
    408             psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    409           }
     400          for (int j = 0; j <= fitOrder; j++) {
     401            if (i + j < 2)
     402              continue;
     403            if (i + j > fitOrder)
     404              continue;
     405            sprintf (name, "PCA1X%1dY%1d", i, j);
     406            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
     407            sprintf (name, "PCA2X%1dY%1d", i, j);
     408            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
     409          }
    410410        }
    411411        psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
    412412      }
    413      
     413
    414414      // remove any existing 'CDi_j style' wcs keywords
    415415      if (psMetadataLookup(header, "CD1_1")) {
    … …  
    419419        psMetadataRemoveKey(header, "CD2_2");
    420420      }
     421
     422      // Remove 'CDi_jX' WCS keywords
     423      psString cd11 = psStringCopy("CD1_1 ");
     424      psString cd12 = psStringCopy("CD1_2 ");
     425      psString cd21 = psStringCopy("CD2_1 ");
     426      psString cd22 = psStringCopy("CD2_2 ");
     427      for (char extra = 'A'; extra <= 'Z'; extra++) {
     428          cd11[strlen(cd11)-1] = extra;
     429          if (psMetadataLookup(header, cd11)) {
     430              cd12[strlen(cd12)-1] = extra;
     431              cd21[strlen(cd21)-1] = extra;
     432              cd22[strlen(cd22)-1] = extra;
     433              psMetadataRemoveKey(header, cd11);
     434              psMetadataRemoveKey(header, cd12);
     435              psMetadataRemoveKey(header, cd21);
     436              psMetadataRemoveKey(header, cd22);
     437          }
     438      }
     439      psFree(cd11);
     440      psFree(cd12);
     441      psFree(cd21);
     442      psFree(cd22);
     443
     444
    421445    } else {
    422446
    … …  
    427451
    428452      if (psMetadataLookup(header, "PC001001")) {
    429         psMetadataRemoveKey(header, "PC001001");
    430         psMetadataRemoveKey(header, "PC001002");
    431         psMetadataRemoveKey(header, "PC002001");
    432         psMetadataRemoveKey(header, "PC002002");
     453        psMetadataRemoveKey(header, "PC001001");
     454        psMetadataRemoveKey(header, "PC001002");
     455        psMetadataRemoveKey(header, "PC002001");
     456        psMetadataRemoveKey(header, "PC002002");
    433457      }
    434458    }
    … …  
    685709    }
    686710
    687     // generate transform with the original orientation (does this rotate about 'center'?) 
     711    // generate transform with the original orientation (does this rotate about 'center'?)
    688712    psPlaneTransform *tpa2 = psPlaneTransformRotate (NULL, tpa1, -1.0*angle);
    689713
    … …  
    967991
    968992    // outFPA projection must be defined as the goal
    969    
     993
    970994    // the output transformations are:
    971995    // chip -> FPA : standard linear trans with needed rotation, etc
    … …  
    9901014        for (int i =  0; i < nSamples; i++) {
    9911015
    992             psSphere srcSky;
    993             psPlane *srcChip = psPlaneAlloc();
    994             psPlane *dstTP = psPlaneAlloc();
     1016            psSphere srcSky;
     1017            psPlane *srcChip = psPlaneAlloc();
     1018            psPlane *dstTP = psPlaneAlloc();
    9951019
    9961020            srcChip->x = bounds->x0 + (i * deltaX / nSamples);
    9971021            srcChip->y = y;
    9981022
    999             psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
    1000             psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
    1001             psDeproject (&srcSky, &srcTP, inFPA->toSky);
    1002            
    1003             // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
    1004 
    1005             psProject (dstTP, &srcSky, outFPA->toSky);
     1023            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
     1024            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
     1025            psDeproject (&srcSky, &srcTP, inFPA->toSky);
     1026
     1027            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
     1028
     1029            psProject (dstTP, &srcSky, outFPA->toSky);
    10061030
    10071031            srcChip->x -= bounds->x0;
    10081032            srcChip->y -= bounds->y0;
    1009             psArrayAdd (src, 100, srcChip);
    1010             psArrayAdd (dst, 100, dstTP);
     1033            psArrayAdd (src, 100, srcChip);
     1034            psArrayAdd (dst, 100, dstTP);
    10111035
    10121036            psFree(srcChip);  // drop our refs to s and d
    … …  
    10211045    if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) {
    10221046        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
    1023         psFree(src);
    1024         psFree(dst);
     1047        psFree(src);
     1048        psFree(dst);
    10251049        return NULL;
    10261050    }
    1027    
     1051
    10281052# if (0)
    10291053    for (int i = 0; i < src->n; i++) {
    1030        
    1031         psSphere srcSky, dstSky;
    1032         psPlane *srcChip = src->data[i];
    1033         psPlane *dstTP   = dst->data[i];
    1034 
    1035         psPlaneTransformApply (&srcFP, newToFPA, srcChip);
    1036         psDeproject (&srcSky, &srcFP, outFPA->toSky);
    1037         psDeproject (&dstSky, dstTP, outFPA->toSky);
    1038 
    1039         double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
    1040         double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
    1041         fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
     1054
     1055        psSphere srcSky, dstSky;
     1056        psPlane *srcChip = src->data[i];
     1057        psPlane *dstTP   = dst->data[i];
     1058
     1059        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
     1060        psDeproject (&srcSky, &srcFP, outFPA->toSky);
     1061        psDeproject (&dstSky, dstTP, outFPA->toSky);
     1062
     1063        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
     1064        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
     1065        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
    10421066
    10431067    }

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

    r28111 r28484  
    759759    psString realName = pmConfigConvertFilename(file->filename, config, create, false);
    760760    if (!realName) {
    761         psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename);
     761        psError(psErrorCodeLast(), false, "failed to determine real name from template for %s\n",
     762                file->filename);
    762763        return false;
    763764    }

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

    r26651 r28484  
    303303
    304304                if (!psThreadJobAddPending(job)) {
    305                     psFree(job);
    306305                    psFree(groups);
    307306                    return false;
    308307                }
    309                 psFree(job);
    310308            }
    311309            if (!psThreadPoolWait(true)) {

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

    r27161 r28484  
    306306    *config = psMetadataConfigRead(NULL, &numBadLines, realName, true);
    307307    if (numBadLines > 0) {
    308         psError(PM_ERR_CONFIG, false, "%d bad lines in %s configuration file (%s)",
     308        psError(PM_ERR_CONFIG, true, "%d bad lines in %s configuration file (%s)",
    309309                numBadLines, description, realName);
    310310        psFree (realName);
    … …  
    17511751    psMetadataItem *item = psMetadataLookup(camera, "FILERULES"); // Item with the file rule of interest
    17521752    if (!item) {
    1753         psError(PM_ERR_CONFIG, false, "Unable to find FILERULES in the camera configuration.");
     1753        psError(PM_ERR_CONFIG, true, "Unable to find FILERULES in the camera configuration.");
    17541754        return NULL;
    17551755    }
    … …  
    17711771    }
    17721772
    1773     return psMetadataLookupMetadata(&mdok, filerules, realname);
     1773    return psMetadataLookupMetadata(NULL, filerules, realname);
    17741774}
    17751775

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

    r25828 r28484  
    88
    99#include "pmConfigMask.h"
     10
     11// Structure to hold the properties of a mask value
     12typedef struct {
     13    char *badMaskName;                  // name for "bad" (i.e., mask me please) pixels
     14    char *fallbackName;                 // Fallback name in case a bad mask name is not defined
     15    psImageMaskType defaultMaskValue;   // Default value in case a bad mask name and its fallback are not defined
     16    bool isBad; // include this value as part of the MASK.VALUE entry (generically bad)
     17} pmConfigMaskInfo;
    1018
    1119static pmConfigMaskInfo masks[] = {

  • branches/pap/psModules/src/config/pmConfigMask.h

    r21183 r28484  
    2020/// @{
    2121
    22 // structure to hold the properties of a mask value
    23 typedef struct {
    24     char *badMaskName;                  // name for "bad" (i.e., mask me please) pixels
    25     char *fallbackName;                 // Fallback name in case a bad mask name is not defined
    26     psImageMaskType defaultMaskValue;   // Default value in case a bad mask name and its fallback are not defined
    27     bool isBad; // include this value as part of the MASK.VALUE entry (generically bad)
    28 } pmConfigMaskInfo;
    29 
    3022// pmConfigMaskSetInMetadata examines named mask values and set the bits for maskValue and
    3123// markValue.  Ensures that the below-named mask values are set, and calculates the mask value
    … …  
    3325// name is not found, or the default values if the fallback name is not found.
    3426bool pmConfigMaskSetInMetadata(psImageMaskType *outMaskValue, // Value of MASK.VALUE, returned
    35                                psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned
    36                                psMetadata *source  // Source of mask bits
     27                               psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned
     28                               psMetadata *source  // Source of mask bits
    3729  );
    3830
    … …  
    4032// Get a mask value by name(s)
    4133psImageMaskType pmConfigMaskGetFromMetadata(psMetadata *source, // Source of masks
    42                                             const char *masks // Mask values to get
     34                                            const char *masks // Mask values to get
    4335  );
    4436
    4537
    46 // lookup an image mask value by name from a psMetadata, without requiring the entry to 
     38// lookup an image mask value by name from a psMetadata, without requiring the entry to
    4739// be of type psImageMaskType, but verifying that it will fit in psImageMaskType
    4840psImageMaskType psMetadataLookupImageMaskFromGeneric (bool *status, const psMetadata *md, const char *name);
    … …  
    5042// Remove from the header keywords starting with the provided string
    5143int pmConfigMaskRemoveHeaderKeywords(psMetadata *header, // Header from which to remove keywords
    52                                      const char *start // Remove keywords that start with this string
     44                                     const char *start // Remove keywords that start with this string
    5345  );
    5446

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

    r21183 r28484  
    144144
    145145            if (!psThreadJobAddPending(job)) {
    146                 psFree(job);
    147146                return false;
    148147            }
    149             psFree(job);
    150148        } else if (!pmBiasSubtractScan(in, sub, scale, xOffset, yOffset, rowStart, rowStop)) {
    151149            psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");

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

    r27597 r28484  
    587587
    588588            if (!psThreadJobAddPending (job)) {
    589                 psFree(job);
    590589                psFree(orders);
    591590                psFree(values);
    592591                return false;
    593592            }
    594             psFree(job);
    595593        } else if (!pmDarkApplyScan(readout, dark, orders, values, bad, doNorm, norm, rowStart, rowStop)) {
    596594            psError(PS_ERR_UNKNOWN, false, "Unable to apply dark.");

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

    r21533 r28484  
    150150
    151151          if (!psThreadJobAddPending(job)) {
    152               psFree(job);
    153152              return false;
    154153          }
    155           psFree(job);
    156154      } else if (!pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat,
    157155                                  xOffset, yOffset, rowStart, rowStop)) {

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

    r27832 r28484  
    351351    psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    352352    if (!psVectorStats (rawStats, counts, NULL, NULL, 0)) {
    353         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    354         return NULL;
     353        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     354        return NULL;
    355355    }
    356356    if (!psVectorStats (resStats, resid, NULL, NULL, 0)) {
    357         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    358         return NULL;
     357        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     358        return NULL;
    359359    }
    360360
    … …  
    794794
    795795                if (!psThreadJobAddPending(job)) {
    796                     psFree(job);
    797796                    return false;
    798797                }
    799                 psFree(job);
    800798            } else if (!pmShutterCorrectionApplyScan(image, mask, var, shutterImage, exptime, blank,
    801799                                                     rowStart, rowStop)) {
    … …  
    10171015
    10181016        if (corr) {
    1019             psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
    1020             if (isfinite(corr->offref) && corr->valid) {
    1021                 psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
    1022                 meanRef += corr->offref;
    1023                 numGood++;
    1024             }
     1017            psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
     1018            if (isfinite(corr->offref) && corr->valid) {
     1019                psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
     1020                meanRef += corr->offref;
     1021                numGood++;
     1022            }
    10251023        } else {
    1026             psTrace("psModules.detrend", 5, "failed Shutter correction fit\n");
    1027         }
     1024            psTrace("psModules.detrend", 5, "failed Shutter correction fit\n");
     1025        }
    10281026
    10291027        psFree(corr);

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

    r27331 r28484  
    3333
    3434
    35 // #define TESTING                         // Enable test output
     35//#define TESTING                         // Enable test output
    3636// #define PEAK_NORM                       // Normalise peaks?
    3737#define PEAK_FLUX 1.0e4                 // Peak flux for each source
    … …  
    235235
    236236            // Get the radius
    237             pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source
     237            pmModel *model = pmModelFromPSFforXY(psf, source->peak->xf, source->peak->yf, PEAK_FLUX); // Model for source
    238238            if (!model || (model->flags & MODEL_MASK)) {
    239239                continue;
    … …  
    321321    numFakes = fakes->n;
    322322
    323     if (numFakes == 0.0) {
     323    if (numFakes == 0) {
    324324        psError(PS_ERR_UNKNOWN, false, "No fake sources are suitable for PSF fitting.");
    325325        psFree(fakes);

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

    r27365 r28484  
    289289                    PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32);
    290290                    if (!psThreadJobAddPending(job)) {
    291                         psFree(job);
    292291                        psFree(source);
    293292                        psFree(target);
    294293                        return NULL;
    295294                    }
    296                     psFree(job);
    297295                } else if (!stackRejectGrow(target, source, kernels, numCols, numRows,
    298296                                            i, xSubMax, j, ySubMax, poorFrac)) {

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

    r28150 r28484  
    12911291
    12921292                if (!psThreadJobAddPending(job)) {
    1293                     psFree(job);
    12941293                    return false;
    12951294                }
    1296                 psFree(job);
    12971295            } else {
    12981296                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,

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

    r27335 r28484  
    860860            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
    861861            if (!psThreadJobAddPending(job)) {
    862                 psFree(job);
    863862                return false;
    864863            }
    865             psFree(job);
    866864        } else {
    867865            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);

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

    r27365 r28484  
    10601060            PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32);
    10611061            if (!psThreadJobAddPending(job)) {
    1062                 psFree(job);
    10631062                return false;
    10641063            }
    1065             psFree(job);
    10661064        } else {
    10671065            if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) {

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

    r27531 r28484  
    174174    PS_ASSERT_PTR_NON_NULL(chip, false);
    175175
    176     if (!pmPSFmodelWrite(chip->analysis, view, file, config)) {
     176    // We need the readout as well, because that has the PSF analysis data (e.g., clumps)
     177    // There is only one, because photometry is done on chip-mosaicked data.
     178    pmFPAview *roView = pmFPAviewAlloc(0); // View to readout
     179    *roView = *view;
     180    roView->cell = 0;
     181    roView->readout = 0;
     182    pmReadout *ro = pmFPAviewThisReadout(roView, chip->parent); // Readout with analysis data
     183    psFree(roView);
     184
     185    if (!pmPSFmodelWrite(chip->analysis, ro ? ro->analysis : NULL, view, file, config)) {
    177186        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
    178187        return false;
    … …  
    189198// else
    190199//   - psf table (+header) : FITS Table
    191 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view,
     200bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view,
    192201                      pmFPAfile *file, pmConfig *config)
    193202{
    … …  
    198207    char *headName, *tableName, *residName;
    199208
    200     if (!analysis) {
     209    if (!chipAnalysis) {
    201210        psError(PM_ERR_PROG, true, "No analysis metadata for chip.");
    202211        return false;
     212    }
     213    if (!roAnalysis) {
     214        psWarning("No analysis metadata for PSF, clump parameters cannot be saved.");
    203215    }
    204216
    … …  
    314326
    315327    // select the psf of interest
    316     pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     328    pmPSF *psf = psMetadataLookupPtr (&status, chipAnalysis, "PSPHOT.PSF");
    317329    if (!psf) {
    318330        psError(PM_ERR_PROG, true, "missing PSF for this analysis metadata");
    … …  
    346358
    347359        // we now save clump parameters for each region : need to save all of those
    348         int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    349         psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions);
    350         for (int i = 0; i < nRegions; i++) {
    351             char regionName[64];
    352             snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    353             psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
    354 
    355             psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   assert (status);
    356             psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   assert (status);
    357             psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  assert (status);
    358             psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  assert (status);
    359 
    360             char key[16];
    361             snprintf (key, 9, "CLX_%03d", i);
    362             psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X);
    363             snprintf (key, 9, "CLY_%03d", i);
    364             psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y);
    365             snprintf (key, 9, "CLDX_%03d", i);
    366             psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX);
    367             snprintf (key, 9, "CLDY_%03d", i);
    368             psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY);
     360        if (roAnalysis) {
     361            int nRegions = psMetadataLookupS32 (&status, roAnalysis, "PSF.CLUMP.NREGIONS");
     362            psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions);
     363            for (int i = 0; i < nRegions; i++) {
     364                char regionName[64];
     365                snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     366                psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
     367
     368                psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   assert (status);
     369                psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   assert (status);
     370                psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  assert (status);
     371                psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  assert (status);
     372
     373                char key[16];
     374                snprintf (key, 9, "CLX_%03d", i);
     375                psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X);
     376                snprintf (key, 9, "CLY_%03d", i);
     377                psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y);
     378                snprintf (key, 9, "CLDX_%03d", i);
     379                psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX);
     380                snprintf (key, 9, "CLDY_%03d", i);
     381                psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY);
     382            }
    369383        }
    370384
    … …  
    492506
    493507        // write the residuals as planes of the image
    494         psArray *images = psArrayAllocEmpty (1);
    495         psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
     508        psArray *images = psArrayAllocEmpty (1);
     509        psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
    496510
    497511        if (psf->residuals->Rx) {
    498512            psArrayAdd (images, 1, psf->residuals->Rx);
    499513            psArrayAdd (images, 1, psf->residuals->Ry);
    500         }
    501 
    502         // note that all N plane are implicitly of the same type, so we convert the mask
    503         if (psf->residuals->mask) {
    504             psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
    505             psArrayAdd (images, 1, mask);
    506             psFree (mask);
    507         }
    508 
    509         // psFitsWriteImageCube (file->fits, header, images, residName);
    510         // psFree (images);
    511 
    512         if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
     514        }
     515
     516        // note that all N plane are implicitly of the same type, so we convert the mask
     517        if (psf->residuals->mask) {
     518            psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     519            psArrayAdd (images, 1, mask);
     520            psFree (mask);
     521        }
     522
     523        // psFitsWriteImageCube (file->fits, header, images, residName);
     524        // psFree (images);
     525
     526        if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
    513527            psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
    514528            psFree(images);
    … …  
    728742    PS_ASSERT_PTR_NON_NULL(file->fpa, false);
    729743
    730     if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
     744    // We need the readout as well, because that has the PSF analysis data (e.g., clumps)
     745    // There may be only one, because photometry is done on chip-mosaicked data.
     746    if (chip->cells->n != 1) {
     747        psError(PM_ERR_PROG, true, "Chip to receive PSF has %ld cells (should be only one)",
     748                chip->cells->n);
     749        return false;
     750    }
     751    pmCell *cell = chip->cells->data[0]; // Cell to receive PSF
     752    pmReadout *ro = NULL;                // Readout to receive PSF
     753    if (cell->readouts->n == 0) {
     754        ro = pmReadoutAlloc(cell);
     755        psFree(ro);                     // Drop reference
     756    } else if (cell->readouts->n != 1) {
     757        psError(PM_ERR_PROG, true, "Cell to receive PSF has %ld readouts (should be only one)",
     758                cell->readouts->n);
     759        return false;
     760    } else {
     761        ro = cell->readouts->data[0];
     762    }
     763    PM_ASSERT_READOUT_NON_NULL(ro, false);
     764    if (!ro->analysis) {
     765        ro->analysis = psMetadataAlloc();
     766    }
     767
     768    if (!pmPSFmodelRead(chip->analysis, ro->analysis, view, file, config)) {
    731769        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
    732770        return false;
    … …  
    737775// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters,
    738776// and header + image for the PSF residual images
    739 bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
     777bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    740778{
     779    PS_ASSERT_METADATA_NON_NULL(chipAnalysis, false);
    741780    PS_ASSERT_PTR_NON_NULL(view, false);
    742781    PS_ASSERT_PTR_NON_NULL(file, false);
    … …  
    803842
    804843    // read the psf clump data for each region
    805     int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN");
    806     if (!status) {
    807         // read old-style psf clump data
    808 
    809         char regionName[64];
    810         snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
    811         psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
    812 
    813         psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
    814         if (!regionMD) {
    815             regionMD = psMetadataAlloc();
    816             psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
    817             psFree (regionMD);
    818         }
    819 
    820         // psf clump data
    821         pmPSFClump psfClump;
    822 
    823         psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  assert(status);
    824         psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  assert(status);
    825         psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  assert(status);
    826         psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  assert(status);
    827 
    828         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    829         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
    830         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    831         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    832     } else {
    833         psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
    834 
    835         for (int i = 0; i < nRegions; i++) {
    836             char key[10];
     844    if (roAnalysis) {
     845        int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN");
     846        if (!status) {
     847            // read old-style psf clump data
     848
    837849            char regionName[64];
    838             snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    839 
    840             psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     850            snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
     851            psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
     852
     853            psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
    841854            if (!regionMD) {
    842855                regionMD = psMetadataAlloc();
    843                 psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
     856                psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
    844857                psFree (regionMD);
    845858            }
    … …  
    848861            pmPSFClump psfClump;
    849862
    850             snprintf (key, 9, "CLX_%03d", i);
    851             psfClump.X  = psMetadataLookupF32 (&status, header, key);  assert(status);
    852             snprintf (key, 9, "CLY_%03d", i);
    853             psfClump.Y  = psMetadataLookupF32 (&status, header, key);  assert(status);
    854             snprintf (key, 9, "CLDX_%03d", i);
    855             psfClump.dX = psMetadataLookupF32 (&status, header, key);  assert(status);
    856             snprintf (key, 9, "CLDY_%03d", i);
    857             psfClump.dY = psMetadataLookupF32 (&status, header, key);  assert(status);
     863            psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  assert(status);
     864            psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  assert(status);
     865            psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  assert(status);
     866            psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  assert(status);
    858867
    859868            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    … …  
    861870            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    862871            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     872        } else {
     873            psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
     874
     875            for (int i = 0; i < nRegions; i++) {
     876                char key[10];
     877                char regionName[64];
     878                snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     879
     880                psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName);
     881                if (!regionMD) {
     882                    regionMD = psMetadataAlloc();
     883                    psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
     884                    psFree (regionMD);
     885                }
     886
     887                // psf clump data
     888                pmPSFClump psfClump;
     889
     890                snprintf (key, 9, "CLX_%03d", i);
     891                psfClump.X  = psMetadataLookupF32 (&status, header, key);  assert(status);
     892                snprintf (key, 9, "CLY_%03d", i);
     893                psfClump.Y  = psMetadataLookupF32 (&status, header, key);  assert(status);
     894                snprintf (key, 9, "CLDX_%03d", i);
     895                psfClump.dX = psMetadataLookupF32 (&status, header, key);  assert(status);
     896                snprintf (key, 9, "CLDY_%03d", i);
     897                psfClump.dY = psMetadataLookupF32 (&status, header, key);  assert(status);
     898
     899                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
     900                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
     901                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
     902                psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     903            }
    863904        }
    864905    }
    … …  
    10191060        }
    10201061
    1021         // note that all N plane are implicitly of the same type, so we convert the mask
    1022         psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
    1023         psImageInit (psf->residuals->mask, 0);
    1024         psImageInit (psf->residuals->Rx, 0.0);
    1025         psImageInit (psf->residuals->Ry, 0.0);
    1026         switch (Nz) {
    1027           case 1: // Ro only
    1028             break;
    1029           case 2: // Ro and mask
     1062        // note that all N plane are implicitly of the same type, so we convert the mask
     1063        psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
     1064        psImageInit (psf->residuals->mask, 0);
     1065        psImageInit (psf->residuals->Rx, 0.0);
     1066        psImageInit (psf->residuals->Ry, 0.0);
     1067        switch (Nz) {
     1068          case 1: // Ro only
     1069            break;
     1070          case 2: // Ro and mask
    10301071            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) {
    10311072                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    10321073                return false;
    10331074            }
    1034             psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
    1035             break;
    1036           case 3: // Ro, Rx and Ry, no mask
     1075            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1076            break;
     1077          case 3: // Ro, Rx and Ry, no mask
    10371078            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
    10381079                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    … …  
    10431084                return false;
    10441085            }
    1045             break;
    1046           case 4: // Ro, Rx, Ry, and mask:
     1086            break;
     1087          case 4: // Ro, Rx, Ry, and mask:
    10471088            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
    10481089                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
    … …  
    10571098                return false;
    10581099            }
    1059             psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
    1060             break;
    1061         }
    1062         psFree (mask);
    1063     }
    1064 
    1065     psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     1100            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
     1101            break;
     1102        }
     1103        psFree (mask);
     1104    }
     1105
     1106    psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
    10661107    psFree (psf);
    10671108

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

    r18601 r28484  
    2020bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    2121bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    22 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
     22bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    2323
    2424bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    … …  
    2727bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    2828bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    29 bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
     29bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
    3030
    3131bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file);

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

    r27531 r28484  
    107107        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    108108        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
    109         psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
    110         psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
    111         source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
    112         source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
    113         source->psfMag = -2.5*log10(source->psfFlux);
     109        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
     110        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
     111        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
     112        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
     113        source->psfMag = -2.5*log10(source->psfFlux);
    114114    } else {
    115115        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
    116         source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     116        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
    117117    }
    118118
    … …  
    120120    if (source->modelFits) {
    121121        bool foundEXT = false;
    122         for (int i = 0; i < source->modelFits->n; i++) {
    123             pmModel *model = source->modelFits->data[i];
    124             status = pmSourcePhotometryModel (&model->mag, NULL, model);
    125             if (model == source->modelEXT) foundEXT = true;
    126         }
    127         if (foundEXT) {
    128             source->extMag = source->modelEXT->mag;
    129         } else {
    130             status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    131         }
     122        for (int i = 0; i < source->modelFits->n; i++) {
     123            pmModel *model = source->modelFits->data[i];
     124            status = pmSourcePhotometryModel (&model->mag, NULL, model);
     125            if (model == source->modelEXT) foundEXT = true;
     126        }
     127        if (foundEXT) {
     128            source->extMag = source->modelEXT->mag;
     129        } else {
     130            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     131        }
    132132    } else {
    133         if (source->modelEXT) {
    134             status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
    135         }
     133        if (source->modelEXT) {
     134            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     135        }
    136136    }
    137137
    … …  
    204204        }
    205205        if (mode & PM_SOURCE_PHOT_APCORR) {
    206             // XXX this should be removed -- we no longer fit for the 'sky bias'
     206            // XXX this should be removed -- we no longer fit for the 'sky bias'
    207207            rflux   = pow (10.0, 0.4*source->psfMag);
    208208            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    … …  
    236236    flux = model->modelFlux (model->params);
    237237    if (flux > 0) {
    238         mag = -2.5*log10(flux);
     238        mag = -2.5*log10(flux);
    239239    }
    240240    if (fitMag) {
    241         *fitMag = mag;
     241        *fitMag = mag;
    242242    }
    243243    if (fitFlux) {
    244         *fitFlux = flux;
     244        *fitFlux = flux;
    245245    }
    246246
    … …  
    380380
    381381    if (source->diffStats == NULL) {
    382         source->diffStats = pmSourceDiffStatsAlloc();
     382        source->diffStats = pmSourceDiffStatsAlloc();
    383383    }
    384384
    … …  
    388388    int   nMask = 0;
    389389    int   nBad  = 0;
    390    
     390
    391391    psImage *flux     = source->pixels;
    392392    psImage *variance = source->variance;
    … …  
    394394
    395395    for (int iy = 0; iy < flux->numRows; iy++) {
    396         for (int ix = 0; ix < flux->numCols; ix++) {
     396        for (int ix = 0; ix < flux->numCols; ix++) {
    397397            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
    398                 nMask ++;
    399                 continue;
    400             }
    401 
    402             float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    403 
    404             if (SN > +FLUX_LIMIT) {
    405                 nGood ++;
    406                 fGood += fabs(flux->data.F32[iy][ix]);
    407             }
    408 
    409             if (SN < -FLUX_LIMIT) {
    410                 nBad ++;
    411                 fBad += fabs(flux->data.F32[iy][ix]);
    412             }
    413         }
     398                nMask ++;
     399                continue;
     400            }
     401
     402            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     403
     404            if (SN > +FLUX_LIMIT) {
     405                nGood ++;
     406                fGood += fabs(flux->data.F32[iy][ix]);
     407            }
     408
     409            if (SN < -FLUX_LIMIT) {
     410                nBad ++;
     411                fBad += fabs(flux->data.F32[iy][ix]);
     412            }
     413        }
    414414    }
    415415
    416416    source->diffStats->nGood      = nGood;
    417     source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);         
    418     source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);         
    419     source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
     417    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);
     418    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);
     419    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);
    420420    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
    421421
    … …  
    628628        }
    629629    }
    630     model->nPix = Npix; 
     630    model->nPix = Npix;
    631631    model->nDOF = Npix - 1;
    632632    model->chisq = dC;
    … …  
    636636
    637637
    638 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
     638double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    639639{
    640640    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    652652    for (int yi = 0; yi < Pi->numRows; yi++) {
    653653        for (int xi = 0; xi < Pi->numCols; xi++) {
    654             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
     654            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    655655                continue;
    656656            if (!unweighted_sum) {
    … …  
    684684}
    685685
    686 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     686double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    687687{
    688688    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    727727    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    728728        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    729             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    730                 continue;
    731             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     729            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     730                continue;
     731            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    732732                continue;
    733733
    … …  
    746746}
    747747
    748 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     748double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    749749{
    750750    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    789789    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    790790        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    791             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    792                 continue;
    793             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     791            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     792                continue;
     793            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    794794                continue;
    795795

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

    r27531 r28484  
    4848    psImage *image,                     ///< image pixels to be used
    4949    psImage *mask,                      ///< mask of pixels to ignore
    50     psImageMaskType maskVal             ///< Value to mask
     50    psImageMaskType maskVal             ///< Value to mask
    5151);
    5252
    … …  
    5858bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
    5959
    60 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    61 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    62 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor);
     60double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     61double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     62double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    6363
    6464// retire these:
Note: See TracChangeset for help on using the changeset viewer.