Index: branches/pap/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 28179) +++ branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 28484) @@ -378,5 +378,5 @@ // XXX make it optional to write out CDi_j terms, or other versions // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity - if (!wcs->wcsCDkeys) { + if (!wcs->wcsCDkeys) { double cdelt1 = wcs->cdelt1; @@ -384,5 +384,5 @@ psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); - + // test the PC00i00j varient: psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1 @@ -390,5 +390,5 @@ psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1 psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2 - + // Elixir-style polynomial terms // XXX currently, Elixir/DVO cannot accept mixed orders @@ -398,18 +398,18 @@ if (fitOrder > 1) { for (int i = 0; i <= fitOrder; i++) { - for (int j = 0; j <= fitOrder; j++) { - if (i + j < 2) - continue; - if (i + j > fitOrder) - continue; - sprintf (name, "PCA1X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - sprintf (name, "PCA2X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - } + for (int j = 0; j <= fitOrder; j++) { + if (i + j < 2) + continue; + if (i + j > fitOrder) + continue; + sprintf (name, "PCA1X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + sprintf (name, "PCA2X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + } } psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder); } - + // remove any existing 'CDi_j style' wcs keywords if (psMetadataLookup(header, "CD1_1")) { @@ -419,4 +419,28 @@ psMetadataRemoveKey(header, "CD2_2"); } + + // Remove 'CDi_jX' WCS keywords + psString cd11 = psStringCopy("CD1_1 "); + psString cd12 = psStringCopy("CD1_2 "); + psString cd21 = psStringCopy("CD2_1 "); + psString cd22 = psStringCopy("CD2_2 "); + for (char extra = 'A'; extra <= 'Z'; extra++) { + cd11[strlen(cd11)-1] = extra; + if (psMetadataLookup(header, cd11)) { + cd12[strlen(cd12)-1] = extra; + cd21[strlen(cd21)-1] = extra; + cd22[strlen(cd22)-1] = extra; + psMetadataRemoveKey(header, cd11); + psMetadataRemoveKey(header, cd12); + psMetadataRemoveKey(header, cd21); + psMetadataRemoveKey(header, cd22); + } + } + psFree(cd11); + psFree(cd12); + psFree(cd21); + psFree(cd22); + + } else { @@ -427,8 +451,8 @@ if (psMetadataLookup(header, "PC001001")) { - psMetadataRemoveKey(header, "PC001001"); - psMetadataRemoveKey(header, "PC001002"); - psMetadataRemoveKey(header, "PC002001"); - psMetadataRemoveKey(header, "PC002002"); + psMetadataRemoveKey(header, "PC001001"); + psMetadataRemoveKey(header, "PC001002"); + psMetadataRemoveKey(header, "PC002001"); + psMetadataRemoveKey(header, "PC002002"); } } @@ -685,5 +709,5 @@ } - // generate transform with the original orientation (does this rotate about 'center'?) + // generate transform with the original orientation (does this rotate about 'center'?) psPlaneTransform *tpa2 = psPlaneTransformRotate (NULL, tpa1, -1.0*angle); @@ -967,5 +991,5 @@ // outFPA projection must be defined as the goal - + // the output transformations are: // chip -> FPA : standard linear trans with needed rotation, etc @@ -990,23 +1014,23 @@ for (int i = 0; i < nSamples; i++) { - psSphere srcSky; - psPlane *srcChip = psPlaneAlloc(); - psPlane *dstTP = psPlaneAlloc(); + psSphere srcSky; + psPlane *srcChip = psPlaneAlloc(); + psPlane *dstTP = psPlaneAlloc(); srcChip->x = bounds->x0 + (i * deltaX / nSamples); srcChip->y = y; - psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip); - psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP); - psDeproject (&srcSky, &srcTP, inFPA->toSky); - - // 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); - - psProject (dstTP, &srcSky, outFPA->toSky); + psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip); + psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP); + psDeproject (&srcSky, &srcTP, inFPA->toSky); + + // 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); + + psProject (dstTP, &srcSky, outFPA->toSky); srcChip->x -= bounds->x0; srcChip->y -= bounds->y0; - psArrayAdd (src, 100, srcChip); - psArrayAdd (dst, 100, dstTP); + psArrayAdd (src, 100, srcChip); + psArrayAdd (dst, 100, dstTP); psFree(srcChip); // drop our refs to s and d @@ -1021,23 +1045,23 @@ if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) { psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed"); - psFree(src); - psFree(dst); + psFree(src); + psFree(dst); return NULL; } - + # if (0) for (int i = 0; i < src->n; i++) { - - psSphere srcSky, dstSky; - psPlane *srcChip = src->data[i]; - psPlane *dstTP = dst->data[i]; - - psPlaneTransformApply (&srcFP, newToFPA, srcChip); - psDeproject (&srcSky, &srcFP, outFPA->toSky); - psDeproject (&dstSky, dstTP, outFPA->toSky); - - double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0; - double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0; - 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); + + psSphere srcSky, dstSky; + psPlane *srcChip = src->data[i]; + psPlane *dstTP = dst->data[i]; + + psPlaneTransformApply (&srcFP, newToFPA, srcChip); + psDeproject (&srcSky, &srcFP, outFPA->toSky); + psDeproject (&dstSky, dstTP, outFPA->toSky); + + double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0; + double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0; + 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); } Index: branches/pap/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 28179) +++ branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 28484) @@ -759,5 +759,6 @@ psString realName = pmConfigConvertFilename(file->filename, config, create, false); if (!realName) { - psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename); + psError(psErrorCodeLast(), false, "failed to determine real name from template for %s\n", + file->filename); return false; } Index: branches/pap/psModules/src/camera/pmReadoutFake.c =================================================================== --- branches/pap/psModules/src/camera/pmReadoutFake.c (revision 28179) +++ branches/pap/psModules/src/camera/pmReadoutFake.c (revision 28484) @@ -303,9 +303,7 @@ if (!psThreadJobAddPending(job)) { - psFree(job); psFree(groups); return false; } - psFree(job); } if (!psThreadPoolWait(true)) { Index: branches/pap/psModules/src/config/pmConfig.c =================================================================== --- branches/pap/psModules/src/config/pmConfig.c (revision 28179) +++ branches/pap/psModules/src/config/pmConfig.c (revision 28484) @@ -306,5 +306,5 @@ *config = psMetadataConfigRead(NULL, &numBadLines, realName, true); if (numBadLines > 0) { - psError(PM_ERR_CONFIG, false, "%d bad lines in %s configuration file (%s)", + psError(PM_ERR_CONFIG, true, "%d bad lines in %s configuration file (%s)", numBadLines, description, realName); psFree (realName); @@ -1751,5 +1751,5 @@ psMetadataItem *item = psMetadataLookup(camera, "FILERULES"); // Item with the file rule of interest if (!item) { - psError(PM_ERR_CONFIG, false, "Unable to find FILERULES in the camera configuration."); + psError(PM_ERR_CONFIG, true, "Unable to find FILERULES in the camera configuration."); return NULL; } @@ -1771,5 +1771,5 @@ } - return psMetadataLookupMetadata(&mdok, filerules, realname); + return psMetadataLookupMetadata(NULL, filerules, realname); } Index: branches/pap/psModules/src/config/pmConfigMask.c =================================================================== --- branches/pap/psModules/src/config/pmConfigMask.c (revision 28179) +++ branches/pap/psModules/src/config/pmConfigMask.c (revision 28484) @@ -8,4 +8,12 @@ #include "pmConfigMask.h" + +// Structure to hold the properties of a mask value +typedef struct { + char *badMaskName; // name for "bad" (i.e., mask me please) pixels + char *fallbackName; // Fallback name in case a bad mask name is not defined + psImageMaskType defaultMaskValue; // Default value in case a bad mask name and its fallback are not defined + bool isBad; // include this value as part of the MASK.VALUE entry (generically bad) +} pmConfigMaskInfo; static pmConfigMaskInfo masks[] = { Index: branches/pap/psModules/src/config/pmConfigMask.h =================================================================== --- branches/pap/psModules/src/config/pmConfigMask.h (revision 28179) +++ branches/pap/psModules/src/config/pmConfigMask.h (revision 28484) @@ -20,12 +20,4 @@ /// @{ -// structure to hold the properties of a mask value -typedef struct { - char *badMaskName; // name for "bad" (i.e., mask me please) pixels - char *fallbackName; // Fallback name in case a bad mask name is not defined - psImageMaskType defaultMaskValue; // Default value in case a bad mask name and its fallback are not defined - bool isBad; // include this value as part of the MASK.VALUE entry (generically bad) -} pmConfigMaskInfo; - // pmConfigMaskSetInMetadata examines named mask values and set the bits for maskValue and // markValue. Ensures that the below-named mask values are set, and calculates the mask value @@ -33,6 +25,6 @@ // name is not found, or the default values if the fallback name is not found. bool pmConfigMaskSetInMetadata(psImageMaskType *outMaskValue, // Value of MASK.VALUE, returned - psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned - psMetadata *source // Source of mask bits + psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned + psMetadata *source // Source of mask bits ); @@ -40,9 +32,9 @@ // Get a mask value by name(s) psImageMaskType pmConfigMaskGetFromMetadata(psMetadata *source, // Source of masks - const char *masks // Mask values to get + const char *masks // Mask values to get ); -// lookup an image mask value by name from a psMetadata, without requiring the entry to +// lookup an image mask value by name from a psMetadata, without requiring the entry to // be of type psImageMaskType, but verifying that it will fit in psImageMaskType psImageMaskType psMetadataLookupImageMaskFromGeneric (bool *status, const psMetadata *md, const char *name); @@ -50,5 +42,5 @@ // Remove from the header keywords starting with the provided string int pmConfigMaskRemoveHeaderKeywords(psMetadata *header, // Header from which to remove keywords - const char *start // Remove keywords that start with this string + const char *start // Remove keywords that start with this string ); Index: branches/pap/psModules/src/detrend/pmBias.c =================================================================== --- branches/pap/psModules/src/detrend/pmBias.c (revision 28179) +++ branches/pap/psModules/src/detrend/pmBias.c (revision 28484) @@ -144,8 +144,6 @@ if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else if (!pmBiasSubtractScan(in, sub, scale, xOffset, yOffset, rowStart, rowStop)) { psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction."); Index: branches/pap/psModules/src/detrend/pmDark.c =================================================================== --- branches/pap/psModules/src/detrend/pmDark.c (revision 28179) +++ branches/pap/psModules/src/detrend/pmDark.c (revision 28484) @@ -587,10 +587,8 @@ if (!psThreadJobAddPending (job)) { - psFree(job); psFree(orders); psFree(values); return false; } - psFree(job); } else if (!pmDarkApplyScan(readout, dark, orders, values, bad, doNorm, norm, rowStart, rowStop)) { psError(PS_ERR_UNKNOWN, false, "Unable to apply dark."); Index: branches/pap/psModules/src/detrend/pmFlatField.c =================================================================== --- branches/pap/psModules/src/detrend/pmFlatField.c (revision 28179) +++ branches/pap/psModules/src/detrend/pmFlatField.c (revision 28484) @@ -150,8 +150,6 @@ if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else if (!pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat, xOffset, yOffset, rowStart, rowStop)) { Index: branches/pap/psModules/src/detrend/pmGainTweak.c =================================================================== --- branches/pap/psModules/src/detrend/pmGainTweak.c (revision 28179) +++ (revision ) @@ -1,59 +1,0 @@ -#include -#include - -#include "pmFPA.h" - -#include "pmGainTweak.h" - -bool pmGainTweak(const pmCell *cell) -{ - PS_ASSERT_PTR_NON_NULL(cell, false); - pmChip *chip = cell->parent; - PS_ASSERT_PTR_NON_NULL(chip, false); - - int numCells = chip->cells; // Number of cells - psVector *gains = psVectorAlloc(numCells, PS_TYPE_F32); // Gain for each cell - psVector *mask = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK); // Mask for gains - psVectorInit(mask, 0); - int numGood = 0; // Number of good gains - - for (int i = 0; i < numCells; i++) { - pmCell *otherCell = chip->cells->data[i]; // A different cell - if (otherCell == cell) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF; - continue; - } - float gain = psMetadataLookupF32(NULL, otherCell->concepts, "CELL.GAIN"); // Gain for cell - gains->data.F32[i] = gain; - if (!isfinite(gains->data.F32[i])) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF; - } else { - numGood++; - } - } - - if (numGood == 0) { - psError(PM_ERR_DETREND, true, "No other cell gains available with which to tweak gain."); - psFree(gains); - psFree(mask); - return false; - } - - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics - if (!psVectorStats(stats, gains, NULL, mask, 0xFF)) { - psError(psErrorCodeLast(), false, "Unable to determine target gain"); - psFree(stats); - psFree(gains); - psFree(mask); - return false; - } - float target = stats->sampleMedian; // Target gain - psFree(stats); - psFree(gains); - psFree(mask); - - psMetadataItem *item = psMetadataLookup(cell->concepts, "CELL.GAIN"); - item->data.F32 = target; - - return true; -} Index: branches/pap/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 28179) +++ branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 28484) @@ -351,10 +351,10 @@ psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); if (!psVectorStats (rawStats, counts, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return NULL; + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return NULL; } if (!psVectorStats (resStats, resid, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return NULL; + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return NULL; } @@ -794,8 +794,6 @@ if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else if (!pmShutterCorrectionApplyScan(image, mask, var, shutterImage, exptime, blank, rowStart, rowStop)) { @@ -1017,13 +1015,13 @@ if (corr) { - psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); - if (isfinite(corr->offref) && corr->valid) { - psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); - meanRef += corr->offref; - numGood++; - } + psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); + if (isfinite(corr->offref) && corr->valid) { + psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); + meanRef += corr->offref; + numGood++; + } } else { - psTrace("psModules.detrend", 5, "failed Shutter correction fit\n"); - } + psTrace("psModules.detrend", 5, "failed Shutter correction fit\n"); + } psFree(corr); Index: branches/pap/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/pap/psModules/src/imcombine/pmPSFEnvelope.c (revision 28179) +++ branches/pap/psModules/src/imcombine/pmPSFEnvelope.c (revision 28484) @@ -33,5 +33,5 @@ -// #define TESTING // Enable test output +//#define TESTING // Enable test output // #define PEAK_NORM // Normalise peaks? #define PEAK_FLUX 1.0e4 // Peak flux for each source @@ -235,5 +235,5 @@ // Get the radius - pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source + pmModel *model = pmModelFromPSFforXY(psf, source->peak->xf, source->peak->yf, PEAK_FLUX); // Model for source if (!model || (model->flags & MODEL_MASK)) { continue; @@ -321,5 +321,5 @@ numFakes = fakes->n; - if (numFakes == 0.0) { + if (numFakes == 0) { psError(PS_ERR_UNKNOWN, false, "No fake sources are suitable for PSF fitting."); psFree(fakes); Index: branches/pap/psModules/src/imcombine/pmStackReject.c =================================================================== --- branches/pap/psModules/src/imcombine/pmStackReject.c (revision 28179) +++ branches/pap/psModules/src/imcombine/pmStackReject.c (revision 28484) @@ -289,10 +289,8 @@ PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32); if (!psThreadJobAddPending(job)) { - psFree(job); psFree(source); psFree(target); return NULL; } - psFree(job); } else if (!stackRejectGrow(target, source, kernels, numCols, numRows, i, xSubMax, j, ySubMax, poorFrac)) { Index: branches/pap/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 28179) +++ branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 28484) @@ -1291,8 +1291,6 @@ if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else { subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, Index: branches/pap/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 28179) +++ branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 28484) @@ -860,8 +860,6 @@ PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else { pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode); Index: branches/pap/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionMatch.c (revision 28179) +++ branches/pap/psModules/src/imcombine/pmSubtractionMatch.c (revision 28484) @@ -1060,8 +1060,6 @@ PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32); if (!psThreadJobAddPending(job)) { - psFree(job); return false; } - psFree(job); } else { if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) { Index: branches/pap/psModules/src/objects/pmPSF_IO.c =================================================================== --- branches/pap/psModules/src/objects/pmPSF_IO.c (revision 28179) +++ branches/pap/psModules/src/objects/pmPSF_IO.c (revision 28484) @@ -174,5 +174,14 @@ PS_ASSERT_PTR_NON_NULL(chip, false); - if (!pmPSFmodelWrite(chip->analysis, view, file, config)) { + // We need the readout as well, because that has the PSF analysis data (e.g., clumps) + // There is only one, because photometry is done on chip-mosaicked data. + pmFPAview *roView = pmFPAviewAlloc(0); // View to readout + *roView = *view; + roView->cell = 0; + roView->readout = 0; + pmReadout *ro = pmFPAviewThisReadout(roView, chip->parent); // Readout with analysis data + psFree(roView); + + if (!pmPSFmodelWrite(chip->analysis, ro ? ro->analysis : NULL, view, file, config)) { psError(psErrorCodeLast(), false, "Failed to write PSF for chip"); return false; @@ -189,5 +198,5 @@ // else // - psf table (+header) : FITS Table -bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, +bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config) { @@ -198,7 +207,10 @@ char *headName, *tableName, *residName; - if (!analysis) { + if (!chipAnalysis) { psError(PM_ERR_PROG, true, "No analysis metadata for chip."); return false; + } + if (!roAnalysis) { + psWarning("No analysis metadata for PSF, clump parameters cannot be saved."); } @@ -314,5 +326,5 @@ // select the psf of interest - pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF"); + pmPSF *psf = psMetadataLookupPtr (&status, chipAnalysis, "PSPHOT.PSF"); if (!psf) { psError(PM_ERR_PROG, true, "missing PSF for this analysis metadata"); @@ -346,25 +358,27 @@ // we now save clump parameters for each region : need to save all of those - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); - psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions); - for (int i = 0; i < nRegions; i++) { - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - - psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); assert (status); - psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); assert (status); - psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); assert (status); - psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); assert (status); - - char key[16]; - snprintf (key, 9, "CLX_%03d", i); - psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X); - snprintf (key, 9, "CLY_%03d", i); - psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y); - snprintf (key, 9, "CLDX_%03d", i); - psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX); - snprintf (key, 9, "CLDY_%03d", i); - psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY); + if (roAnalysis) { + int nRegions = psMetadataLookupS32 (&status, roAnalysis, "PSF.CLUMP.NREGIONS"); + psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions); + for (int i = 0; i < nRegions; i++) { + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); + + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); assert (status); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); assert (status); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); assert (status); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); assert (status); + + char key[16]; + snprintf (key, 9, "CLX_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X); + snprintf (key, 9, "CLY_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y); + snprintf (key, 9, "CLDX_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX); + snprintf (key, 9, "CLDY_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY); + } } @@ -492,23 +506,23 @@ // write the residuals as planes of the image - psArray *images = psArrayAllocEmpty (1); - psArrayAdd (images, 1, psf->residuals->Ro); // z = 0 is Ro + psArray *images = psArrayAllocEmpty (1); + psArrayAdd (images, 1, psf->residuals->Ro); // z = 0 is Ro if (psf->residuals->Rx) { psArrayAdd (images, 1, psf->residuals->Rx); psArrayAdd (images, 1, psf->residuals->Ry); - } - - // note that all N plane are implicitly of the same type, so we convert the mask - if (psf->residuals->mask) { - psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type); - psArrayAdd (images, 1, mask); - psFree (mask); - } - - // psFitsWriteImageCube (file->fits, header, images, residName); - // psFree (images); - - if (!psFitsWriteImageCube (file->fits, header, images, residName)) { + } + + // note that all N plane are implicitly of the same type, so we convert the mask + if (psf->residuals->mask) { + psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type); + psArrayAdd (images, 1, mask); + psFree (mask); + } + + // psFitsWriteImageCube (file->fits, header, images, residName); + // psFree (images); + + if (!psFitsWriteImageCube (file->fits, header, images, residName)) { psError(psErrorCodeLast(), false, "Unable to write PSF residuals."); psFree(images); @@ -728,5 +742,29 @@ PS_ASSERT_PTR_NON_NULL(file->fpa, false); - if (!pmPSFmodelRead (chip->analysis, view, file, config)) { + // We need the readout as well, because that has the PSF analysis data (e.g., clumps) + // There may be only one, because photometry is done on chip-mosaicked data. + if (chip->cells->n != 1) { + psError(PM_ERR_PROG, true, "Chip to receive PSF has %ld cells (should be only one)", + chip->cells->n); + return false; + } + pmCell *cell = chip->cells->data[0]; // Cell to receive PSF + pmReadout *ro = NULL; // Readout to receive PSF + if (cell->readouts->n == 0) { + ro = pmReadoutAlloc(cell); + psFree(ro); // Drop reference + } else if (cell->readouts->n != 1) { + psError(PM_ERR_PROG, true, "Cell to receive PSF has %ld readouts (should be only one)", + cell->readouts->n); + return false; + } else { + ro = cell->readouts->data[0]; + } + PM_ASSERT_READOUT_NON_NULL(ro, false); + if (!ro->analysis) { + ro->analysis = psMetadataAlloc(); + } + + if (!pmPSFmodelRead(chip->analysis, ro->analysis, view, file, config)) { psError(psErrorCodeLast(), false, "Failed to write PSF for chip"); return false; @@ -737,6 +775,7 @@ // for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, // and header + image for the PSF residual images -bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { + PS_ASSERT_METADATA_NON_NULL(chipAnalysis, false); PS_ASSERT_PTR_NON_NULL(view, false); PS_ASSERT_PTR_NON_NULL(file, false); @@ -803,43 +842,17 @@ // read the psf clump data for each region - int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN"); - if (!status) { - // read old-style psf clump data - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); - psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); - - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - if (!regionMD) { - regionMD = psMetadataAlloc(); - psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); - psFree (regionMD); - } - - // psf clump data - pmPSFClump psfClump; - - psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); - psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); - psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); - psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); - - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - } else { - psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); - - for (int i = 0; i < nRegions; i++) { - char key[10]; + if (roAnalysis) { + int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN"); + if (!status) { + // read old-style psf clump data + char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); + psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); if (!regionMD) { regionMD = psMetadataAlloc(); - psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); psFree (regionMD); } @@ -848,12 +861,8 @@ pmPSFClump psfClump; - snprintf (key, 9, "CLX_%03d", i); - psfClump.X = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLY_%03d", i); - psfClump.Y = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLDX_%03d", i); - psfClump.dX = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLDY_%03d", i); - psfClump.dY = psMetadataLookupF32 (&status, header, key); assert(status); + psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); + psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); + psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); + psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); @@ -861,4 +870,36 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + } else { + psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); + + for (int i = 0; i < nRegions; i++) { + char key[10]; + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + + psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + + // psf clump data + pmPSFClump psfClump; + + snprintf (key, 9, "CLX_%03d", i); + psfClump.X = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLY_%03d", i); + psfClump.Y = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLDX_%03d", i); + psfClump.dX = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLDY_%03d", i); + psfClump.dY = psMetadataLookupF32 (&status, header, key); assert(status); + + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + } } } @@ -1019,20 +1060,20 @@ } - // note that all N plane are implicitly of the same type, so we convert the mask - psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type); - psImageInit (psf->residuals->mask, 0); - psImageInit (psf->residuals->Rx, 0.0); - psImageInit (psf->residuals->Ry, 0.0); - switch (Nz) { - case 1: // Ro only - break; - case 2: // Ro and mask + // note that all N plane are implicitly of the same type, so we convert the mask + psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type); + psImageInit (psf->residuals->mask, 0); + psImageInit (psf->residuals->Rx, 0.0); + psImageInit (psf->residuals->Ry, 0.0); + switch (Nz) { + case 1: // Ro only + break; + case 2: // Ro and mask if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); return false; } - psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); - break; - case 3: // Ro, Rx and Ry, no mask + psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); + break; + case 3: // Ro, Rx and Ry, no mask if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); @@ -1043,6 +1084,6 @@ return false; } - break; - case 4: // Ro, Rx, Ry, and mask: + break; + case 4: // Ro, Rx, Ry, and mask: if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) { psError(psErrorCodeLast(), false, "Unable to read PSF residual image."); @@ -1057,11 +1098,11 @@ return false; } - psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); - break; - } - psFree (mask); - } - - psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK); + break; + } + psFree (mask); + } + + psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); psFree (psf); Index: branches/pap/psModules/src/objects/pmPSF_IO.h =================================================================== --- branches/pap/psModules/src/objects/pmPSF_IO.h (revision 28179) +++ branches/pap/psModules/src/objects/pmPSF_IO.h (revision 28484) @@ -20,5 +20,5 @@ bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, pmConfig *config); bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, pmConfig *config); -bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config); +bool pmPSFmodelWrite (const psMetadata *chipAnalysis, const psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, pmConfig *config); bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); @@ -27,5 +27,5 @@ bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmPSFmodelRead (psMetadata *chipAnalysis, psMetadata *roAnalysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file); Index: branches/pap/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- branches/pap/psModules/src/objects/pmSourcePhotometry.c (revision 28179) +++ branches/pap/psModules/src/objects/pmSourcePhotometry.c (revision 28484) @@ -107,12 +107,12 @@ // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y); - psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); - psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); - source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]; - source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0]; - source->psfMag = -2.5*log10(source->psfFlux); + psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); + psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); + source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]; + source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0]; + source->psfMag = -2.5*log10(source->psfFlux); } else { status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF); - source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]); + source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]); } @@ -120,18 +120,18 @@ if (source->modelFits) { bool foundEXT = false; - for (int i = 0; i < source->modelFits->n; i++) { - pmModel *model = source->modelFits->data[i]; - status = pmSourcePhotometryModel (&model->mag, NULL, model); - if (model == source->modelEXT) foundEXT = true; - } - if (foundEXT) { - source->extMag = source->modelEXT->mag; - } else { - status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); - } + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + status = pmSourcePhotometryModel (&model->mag, NULL, model); + if (model == source->modelEXT) foundEXT = true; + } + if (foundEXT) { + source->extMag = source->modelEXT->mag; + } else { + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); + } } else { - if (source->modelEXT) { - status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); - } + if (source->modelEXT) { + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); + } } @@ -204,5 +204,5 @@ } if (mode & PM_SOURCE_PHOT_APCORR) { - // XXX this should be removed -- we no longer fit for the 'sky bias' + // XXX this should be removed -- we no longer fit for the 'sky bias' rflux = pow (10.0, 0.4*source->psfMag); source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux; @@ -236,11 +236,11 @@ flux = model->modelFlux (model->params); if (flux > 0) { - mag = -2.5*log10(flux); + mag = -2.5*log10(flux); } if (fitMag) { - *fitMag = mag; + *fitMag = mag; } if (fitFlux) { - *fitFlux = flux; + *fitFlux = flux; } @@ -380,5 +380,5 @@ if (source->diffStats == NULL) { - source->diffStats = pmSourceDiffStatsAlloc(); + source->diffStats = pmSourceDiffStatsAlloc(); } @@ -388,5 +388,5 @@ int nMask = 0; int nBad = 0; - + psImage *flux = source->pixels; psImage *variance = source->variance; @@ -394,28 +394,28 @@ for (int iy = 0; iy < flux->numRows; iy++) { - for (int ix = 0; ix < flux->numCols; ix++) { + for (int ix = 0; ix < flux->numCols; ix++) { if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) { - nMask ++; - continue; - } - - float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - - if (SN > +FLUX_LIMIT) { - nGood ++; - fGood += fabs(flux->data.F32[iy][ix]); - } - - if (SN < -FLUX_LIMIT) { - nBad ++; - fBad += fabs(flux->data.F32[iy][ix]); - } - } + nMask ++; + continue; + } + + float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + + if (SN > +FLUX_LIMIT) { + nGood ++; + fGood += fabs(flux->data.F32[iy][ix]); + } + + if (SN < -FLUX_LIMIT) { + nBad ++; + fBad += fabs(flux->data.F32[iy][ix]); + } + } } source->diffStats->nGood = nGood; - source->diffStats->fRatio = (fGood + fBad == 0.0) ? NAN : fGood / (fGood + fBad); - source->diffStats->nRatioBad = (nGood + nBad == 0) ? NAN : nGood / (float) (nGood + nBad); - source->diffStats->nRatioMask = (nGood + nMask == 0) ? NAN : nGood / (float) (nGood + nMask); + source->diffStats->fRatio = (fGood + fBad == 0.0) ? NAN : fGood / (fGood + fBad); + source->diffStats->nRatioBad = (nGood + nBad == 0) ? NAN : nGood / (float) (nGood + nBad); + source->diffStats->nRatioMask = (nGood + nMask == 0) ? NAN : nGood / (float) (nGood + nMask); source->diffStats->nRatioAll = (nGood + nMask + nBad == 0) ? NAN : nGood / (float) (nGood + nMask + nBad); @@ -628,5 +628,5 @@ } } - model->nPix = Npix; + model->nPix = Npix; model->nDOF = Npix - 1; model->chisq = dC; @@ -636,5 +636,5 @@ -double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor) +double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -652,5 +652,5 @@ for (int yi = 0; yi < Pi->numRows; yi++) { for (int xi = 0; xi < Pi->numCols; xi++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) continue; if (!unweighted_sum) { @@ -684,5 +684,5 @@ } -double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor) +double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -727,7 +727,7 @@ for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) - continue; - if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) + continue; + if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal) continue; @@ -746,5 +746,5 @@ } -double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor) +double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(Mi, NAN); @@ -789,7 +789,7 @@ for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi]) - continue; - if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj]) + if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal) + continue; + if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal) continue; Index: branches/pap/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- branches/pap/psModules/src/objects/pmSourcePhotometry.h (revision 28179) +++ branches/pap/psModules/src/objects/pmSourcePhotometry.h (revision 28484) @@ -48,5 +48,5 @@ psImage *image, ///< image pixels to be used psImage *mask, ///< mask of pixels to ignore - psImageMaskType maskVal ///< Value to mask + psImageMaskType maskVal ///< Value to mask ); @@ -58,7 +58,7 @@ bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal); -double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor); -double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor); -double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor); +double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); +double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); +double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); // retire these: