Index: branches/pap/psModules/src/astrom/pmAstrometryDistortion.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryDistortion.c (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryDistortion.c (revision 25027) @@ -151,8 +151,14 @@ // also measure the L and M median positions as a representative coordinate - psVectorStats (stats, L, NULL, NULL, 0); + if (!psVectorStats (stats, L, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + goto skip; + } grad->FP.x = stats->sampleMedian; - psVectorStats (stats, M, NULL, NULL, 0); + if (!psVectorStats (stats, M, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + goto skip; + } grad->FP.y = stats->sampleMedian; Index: branches/pap/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryObjects.c (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryObjects.c (revision 25027) @@ -83,7 +83,12 @@ { psArray *matches = psArrayAllocEmpty(x1->n); + psVector *found1 = psVectorAlloc(x1->n, PS_TYPE_S8); + psVector *found2 = psVectorAlloc(x2->n, PS_TYPE_S8); const double RADIUS_SQR = PS_SQR(RADIUS); double dX, dY, dR; + + psVectorInit (found1, 0); + psVectorInit (found2, 0); int jStart; @@ -100,6 +105,15 @@ } + if (found1->data.S8[i]) { + i++; + continue; + } + if (found2->data.S8[j]) { + j++; + continue; + } + jStart = j; - while (fabs(dX) < RADIUS && j < x2->n) { + while ((fabs(dX) < RADIUS) && (j < x2->n)) { dX = x1->data.F64[i] - x2->data.F64[j]; @@ -111,4 +125,8 @@ continue; } + if (found2->data.S8[j]) { + j++; + continue; + } // got a match; add to output list @@ -117,4 +135,7 @@ psFree (match); + found1->data.S8[i] = 1; + found2->data.S8[j] = 1; + j++; } @@ -122,4 +143,7 @@ i++; } + psFree (found1); + psFree (found2); + return (matches); } @@ -420,4 +444,5 @@ obj->sky = psSphereAlloc(); obj->Mag = 0; + obj->Color= 0; obj->dMag = 0; @@ -447,4 +472,5 @@ *obj->sky = *old->sky; obj->Mag = old->Mag; + obj->Color = old->Color; obj->dMag = old->dMag; Index: branches/pap/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryObjects.h (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryObjects.h (revision 25027) @@ -33,12 +33,13 @@ typedef struct { - psPlane *pix; ///< the position in the pmReadout frame - psPlane *cell; ///< the position in the pmCell frame - psPlane *chip; ///< the position in the pmChip frame - psPlane *FP; ///< the position in the pmFPA frame - psPlane *TP; ///< the position in the tangent plane - psSphere *sky; ///< the position on the Celestial Sphere. - double Mag; ///< object magnitude XXX what filter? - double dMag; ///< error on object magnitude + psPlane *pix; ///< the position in the pmReadout frame + psPlane *cell; ///< the position in the pmCell frame + psPlane *chip; ///< the position in the pmChip frame + psPlane *FP; ///< the position in the pmFPA frame + psPlane *TP; ///< the position in the tangent plane + psSphere *sky; ///< the position on the Celestial Sphere. + float Mag; ///< object magnitude in extracted filter + float Color; ///< object color + float dMag; ///< error on object magnitude } pmAstromObj; Index: branches/pap/psModules/src/astrom/pmAstrometryUtils.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryUtils.c (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryUtils.c (revision 25027) @@ -46,7 +46,7 @@ /* this loop uses the Newton-Raphson method to solve for Xo,Yo - * it needs the high order terms to be small - * Xo,Yo are in pixels; - */ + * it needs the high order terms to be small + * Xo,Yo are in pixels; + */ double dPos = tol + 1; for (int i = 0; (dPos > tol) && (i < 20); i++) { @@ -60,12 +60,12 @@ Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo); - if (!psMatrixGJSolveF32 (Alpha, Beta)) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve for center."); - psFree (Alpha); - psFree (Beta); - psFree (XdX); - psFree (XdY); - psFree (YdX); - psFree (YdY); + if (!psMatrixGJSolve (Alpha, Beta)) { + psError(PS_ERR_UNKNOWN, false, "Unable to solve for center."); + psFree (Alpha); + psFree (Beta); + psFree (XdX); + psFree (XdY); + psFree (YdX); + psFree (YdY); return NULL; } Index: branches/pap/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryVisual.c (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryVisual.c (revision 25027) @@ -17,7 +17,7 @@ #include "pmFPAExtent.h" -# if (HAVE_KAPA) -# include -# include "pmKapaPlots.h" +#if (HAVE_KAPA) +#include +#include "pmKapaPlots.h" #include "pmVisual.h" Index: branches/pap/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 25027) @@ -96,4 +96,9 @@ if (!status1 || !status2) { + Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1"); + Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2"); + } + + if (!status1 || !status2) { Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1"); Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2"); @@ -143,6 +148,6 @@ int Nx = region->x1 - region->x0; int Ny = region->y1 - region->y0; - psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Mosaic Dimensions", Nx); - psMetadataAddS32 (header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Mosaic Dimensions", Ny); + psMetadataAddS32 (header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "Mosaic Dimensions", Nx); + psMetadataAddS32 (header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "Mosaic Dimensions", Ny); pmAstromWCStoHeader (header, wcs); Index: branches/pap/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryWCS.h (revision 23948) +++ branches/pap/psModules/src/astrom/pmAstrometryWCS.h (revision 25027) @@ -68,6 +68,7 @@ bool psPlaneDistortIsDiagonal (psPlaneDistort *distort); -# define PM_DEG_RAD 57.295779513082322 -# define PM_RAD_DEG 0.017453292519943 +// XXX probably should remove these and just use the PS_ version in the code +# define PM_DEG_RAD PS_DEG_RAD +# define PM_RAD_DEG PS_RAD_DEG /// @} Index: branches/pap/psModules/src/camera/Makefile.am =================================================================== --- branches/pap/psModules/src/camera/Makefile.am (revision 23948) +++ branches/pap/psModules/src/camera/Makefile.am (revision 25027) @@ -23,4 +23,5 @@ pmFPAfileIO.c \ pmFPAfileFitsIO.c \ + pmFPAfileFringeIO.c \ pmFPAFlags.c \ pmFPALevel.c \ @@ -51,4 +52,5 @@ pmFPAfileIO.h \ pmFPAfileFitsIO.h \ + pmFPAfileFringeIO.h \ pmFPAFlags.h \ pmFPALevel.h \ Index: branches/pap/psModules/src/camera/pmFPABin.c =================================================================== --- branches/pap/psModules/src/camera/pmFPABin.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPABin.c (revision 25027) @@ -28,10 +28,17 @@ int numColsOut = binning->nXruff, numRowsOut = binning->nYruff; // Size of output image - // Output image - psImage *outImage; + + psImage *outImage; // Output image if (out->image && out->image->numCols >= numColsOut && out->image->numRows >= numRowsOut) { outImage = out->image; } else { outImage = out->image = psImageRecycle(out->image, numColsOut, numRowsOut, PS_TYPE_F32); + } + + psImage *outMask; // Output mask + if (out->mask && out->mask->numCols >= numColsOut && out->mask->numRows >= numRowsOut) { + outMask = out->mask; + } else { + outMask = out->mask = psImageRecycle(out->mask, numColsOut, numRowsOut, PS_TYPE_IMAGE_MASK); } @@ -53,4 +60,7 @@ continue; } + if (!isfinite(inImage->data.F32[y][x])) { + continue; + } sum += inImage->data.F32[y][x]; numPix++; @@ -58,5 +68,16 @@ } - outImage->data.F32[yOut][xOut] = numPix > 0 ? sum / numPix : NAN; + // Values to set + float imageValue; + psImageMaskType maskValue; + if (numPix > 0) { + imageValue = sum / numPix; + maskValue = 0; + } else { + imageValue = NAN; + maskValue = maskVal; + } + outImage->data.F32[yOut][xOut] = imageValue; + outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue; xStart = xStop; } Index: branches/pap/psModules/src/camera/pmFPACopy.c =================================================================== --- branches/pap/psModules/src/camera/pmFPACopy.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPACopy.c (revision 25027) @@ -70,4 +70,5 @@ psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); *target = psImageAlloc(binning->nXruff, binning->nYruff, source->type.type); + psImageInit (*target, 0.0); return; } @@ -221,4 +222,11 @@ psTrace("psModules.camera", 3, "xFlip: %d; yFlip: %d\n", xFlip, yFlip); + // Blow away extant readouts + for (int i = 0; i < target->readouts->n; i++) { + psFree(target->readouts->data[i]); + target->readouts->data[i] = NULL; + } + target->readouts->n = 0; + // Perform deep copy of the images. I would prefer *not* to do a deep copy, in the interests of speed (we // still need to do another deep copy into the HDU for when we write out), but this is the only way I can @@ -242,4 +250,21 @@ readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip, pixels); + // Copy covariance matrix: doesn't care about flips, etc. + if (sourceReadout->covariance) { + if (targetReadout->covariance) { + psFree(targetReadout->covariance); + } + targetReadout->covariance = psKernelCopy(sourceReadout->covariance); +#if 0 + if (binning) { + // XXX This isn't strictly correct, but we don't have a function that bins covariance matrices + // with unequal binning factors. + psKernel *covar = psImageCovarianceBin(PS_MAX(binning->nXbin, binning->nYbin), + targetReadout->covariance); + psFree(targetReadout->covariance); + targetReadout->covariance = covar; + } +#endif + } // Copy bias Index: branches/pap/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAMaskWeight.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPAMaskWeight.c (revision 25027) @@ -109,14 +109,24 @@ float saturation = psMetadataLookupF32(&mdok, cell->concepts, "CELL.SATURATION"); // Saturation level if (!mdok || isnan(saturation)) { - psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n"); - return false; + // psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n"); + // return false; + psWarning("CELL.SATURATION is not set --- completely masking cell.\n"); + saturation = NAN; } float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level if (!mdok || isnan(bad)) { - psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n"); - return false; + // psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n"); + // return false; + psWarning("CELL.BAD is not set --- completely masking cell.\n"); + bad = NAN; } psTrace("psModules.camera", 5, "Saturation: %f, bad: %f\n", saturation, bad); + // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN; + // in this case, we have to set the mask as well + float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain + if (!mdok) { gain = NAN; } + float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise + if (!mdok) { readnoise = NAN; } // Set up the mask @@ -127,4 +137,11 @@ } psImage *mask = readout->mask; // The mask pixels + + // completely mask if SATURATION or BAD are invalid + if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) { + psImageInit(mask, badMask); + return true; + } + psImageInit(mask, 0); @@ -199,7 +216,8 @@ } -bool pmReadoutSetVariance(pmReadout *readout, bool poisson) +bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); + // check that the noiseMap (if it exists) matches the readout variance size) pmCell *cell = readout->parent; // The parent cell @@ -209,15 +227,31 @@ float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain if (!mdok || isnan(gain)) { - psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n"); - return false; + // psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n"); + // return false; + psWarning("CELL.GAIN is not set --- setting variance to NAN\n"); + gain = NAN; } float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise if (!mdok || isnan(readnoise)) { - psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n"); - return false; - } - if (psMetadataLookup(cell->concepts, "CELL.READNOISE.UPDATE")) { + // psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n"); + // return false; + psWarning("CELL.READNOISE is not set --- setting variance to NAN\n"); + readnoise = NAN; + } + // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN) + if (!isnan(gain) && psMetadataLookup(cell->concepts, "CELL.READNOISE.UPDATE")) { psError(PS_ERR_IO, true, "CELL.READNOISE has not yet been updated for the gain"); return false; + } + + // for invalid input data, set the readout variance to NAN + if (isnan(gain) || isnan(readnoise)) { + if (!readout->variance) { + // generate the image if needed + readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); + } + // XXX need to set the mask, if defined + psImageInit(readout->variance, NAN); + return true; } @@ -228,4 +262,5 @@ // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1. + // XXX this calculation is wrong: limit is 1 e-, but this is in DN readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs"); readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max", @@ -239,6 +274,12 @@ } - readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", - psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + // apply a supplied readnoise map (NOTE: in DN, not electrons): + if (noiseMap) { + psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap); + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar); + psFree (rdVar); + } else { + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + } return true; @@ -247,5 +288,5 @@ // this function creates the variance pixels, or uses the existing variance pixels. it will set // the noise pixel values only if the variance image is not supplied -bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson) +bool pmReadoutGenerateVariance(pmReadout *readout, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -291,8 +332,8 @@ readout->variance = variance; - return pmReadoutSetVariance(readout, poisson); -} - -bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson) + return pmReadoutSetVariance(readout, noiseMap, poisson); +} + +bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -301,10 +342,10 @@ success &= pmReadoutGenerateMask(readout, satMask, badMask); - success &= pmReadoutGenerateVariance(readout, poisson); + success &= pmReadoutGenerateVariance(readout, noiseMap, poisson); return success; } -bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson) +bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(cell, false); @@ -314,5 +355,5 @@ for (int i = 0; i < readouts->n; i++) { pmReadout *readout = readouts->data[i]; // The readout - success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask); + success &= pmReadoutGenerateMaskVariance(readout, satMask, badMask, noiseMap, poisson); } @@ -482,5 +523,5 @@ return false; } - float stdev = psStatsGetValue(stdevStats, stdevStat); // Stadard deviation of fluxes + float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes psFree(stdevStats); psFree(noise); Index: branches/pap/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- branches/pap/psModules/src/camera/pmFPAMaskWeight.h (revision 23948) +++ branches/pap/psModules/src/camera/pmFPAMaskWeight.h (revision 25027) @@ -54,4 +54,5 @@ /// can't be generated. bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -72,4 +73,5 @@ /// with HDU entry). This is intended for most operations. bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -81,4 +83,5 @@ psImageMaskType sat, ///< Mask value to give saturated pixels psImageMaskType bad, ///< Mask value to give bad (low) pixels + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -90,4 +93,5 @@ psImageMaskType sat, ///< Mask value to give saturated pixels psImageMaskType bad, ///< Mask value to give bad (low) pixels + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); Index: branches/pap/psModules/src/camera/pmFPARead.c =================================================================== --- branches/pap/psModules/src/camera/pmFPARead.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPARead.c (revision 25027) @@ -64,5 +64,5 @@ // Set the "thisXXXScan" value in the readout for the appropriate image type -static int readoutSetThisScan(pmReadout *readout, // Readout of interest +static void readoutSetThisScan(pmReadout *readout, // Readout of interest fpaReadType type, // Type of image int thisScan // Starting scan number @@ -72,15 +72,14 @@ case FPA_READ_TYPE_IMAGE: readout->thisImageScan = thisScan; - return readout->lastImageScan; + return; case FPA_READ_TYPE_MASK: readout->thisMaskScan = thisScan; - return readout->lastMaskScan; + return; case FPA_READ_TYPE_VARIANCE: readout->thisVarianceScan = thisScan; - return readout->lastVarianceScan; + return; default: psAbort("Unknown read type: %x\n", type); } - return false; } @@ -103,5 +102,5 @@ // Set the "lastXXXScan" value in the readout for the appropriate image type -static int readoutSetLastScan(pmReadout *readout, // Readout of interest +static void readoutSetLastScan(pmReadout *readout, // Readout of interest fpaReadType type, // Type of image int lastScan // Last scan number @@ -111,15 +110,14 @@ case FPA_READ_TYPE_IMAGE: readout->lastImageScan = lastScan; - return readout->lastImageScan; + return; case FPA_READ_TYPE_MASK: readout->lastMaskScan = lastScan; - return readout->lastMaskScan; + return; case FPA_READ_TYPE_VARIANCE: readout->lastVarianceScan = lastScan; - return readout->lastVarianceScan; + return; default: psAbort("Unknown read type: %x\n", type); } - return false; } @@ -143,5 +141,5 @@ // Determine number of readouts in the FITS file // In the process, reads the header and concepts -static bool cellNumReadouts(pmCell *cell, // Cell of interest +static int cellNumReadouts(pmCell *cell, // Cell of interest psFits *fits, // FITS file pmConfig *config // Configuration @@ -155,13 +153,13 @@ if (!hdu || hdu->blankPHU) { psError(PS_ERR_IO, true, "Unable to find HDU"); - return false; + return 0; } if (!pmCellReadHeader(cell, fits, config)) { psError(PS_ERR_IO, false, "Unable to read header for cell!\n"); - return false; + return 0; } if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS, true, NULL)) { psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n"); - return false; + return 0; } @@ -171,15 +169,16 @@ if (!mdok) { psError(PS_ERR_IO, true, "Unable to find NAXIS in header for extension %s\n", hdu->extname); - return false; - } + return 0; + } + if (naxis == 0) { // No pixels to read psError(PS_ERR_IO, true, "No pixels in extension %s.", hdu->extname); - return false; + return 0; } if (naxis < 2 || naxis > 3) { psError(PS_ERR_IO, true, "NAXIS in header of extension %s (= %d) is not valid.\n", hdu->extname, naxis); - return false; + return 0; } int naxis3; // Number of image planes @@ -188,5 +187,5 @@ if (!mdok) { psError(PS_ERR_IO, true, "Unable to find NAXIS3 in header for extension %s\n", hdu->extname); - return false; + return 0; } } else { @@ -296,5 +295,6 @@ static bool readoutMore(pmReadout *readout, // Readout of interest psFits *fits, // FITS file - int z, // Plane number + int z, // Plane number to read + int *zMax, // Max plane number in this cell int numScans, // Number of scans to read at a time fpaReadType type, // Type of image @@ -310,11 +310,10 @@ // N readouts, but numScans set to 0. only the first should report that it requires data, // even if all readouts lack the image pointer. - if (!image) { + if (numScans == 0) { + if (!image) { return true; - } - - // If we have already read an image, this result implies we are done (no more to read) - if (numScans == 0) { - return false; + } else { + return false; + } } @@ -324,6 +323,6 @@ return false; } - int naxis3 = cellNumReadouts(cell, fits, config); // Number of planes - if (z >= naxis3) { + *zMax = cellNumReadouts(cell, fits, config); // Number of planes + if (z >= *zMax) { // No more to read return false; @@ -486,4 +485,5 @@ psFits *fits, // FITS file int z, // Desired image plane + int *zMax, // Max plane number in this cell int numScans, // Number of scans (row or col depends on CELL.READDIR); 0 for all int overlap, // Number of scans (row/col) to overlap between scans @@ -511,4 +511,6 @@ int naxis3 = cellNumReadouts(cell, fits, config); // Number of image planes + if (zMax) *zMax = naxis3; + if (z >= naxis3) { psError(PS_ERR_IO, false, "Desired image plane (%d) exceeds available number (%d).", @@ -592,6 +594,5 @@ } - // Determine the number of scans to read - int lastScan = readoutSetLastScan(readout, type, thisScan + numScans); + int origThisScan = thisScan; // Original value of thisScan (starting point for read) if (thisScan == 0) { overlap = 0; @@ -601,5 +602,4 @@ thisScan = 0; } - readoutSetThisScan(readout, type, thisScan); // Calculate limits, adjust readout->row0,col0 @@ -620,4 +620,8 @@ } } + int lastScan = origThisScan + numScans; // Last scan to read + + readoutSetThisScan(readout, type, thisScan); + readoutSetLastScan(readout, type, lastScan); // Blow away existing data. @@ -1021,13 +1025,13 @@ -bool pmReadoutMore(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) +bool pmReadoutMore(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_IMAGE, config); -} - -bool pmReadoutReadChunk(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, pmConfig *config) + return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_IMAGE, config); +} + +bool pmReadoutReadChunk(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -1036,5 +1040,5 @@ PS_ASSERT_INT_NONNEGATIVE(numScans, false); - return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_IMAGE, config); + return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_IMAGE, config); } @@ -1044,5 +1048,5 @@ PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_IMAGE, config); + return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_IMAGE, config); } @@ -1084,13 +1088,13 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool pmReadoutMoreMask(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) +bool pmReadoutMoreMask(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_MASK, config); -} - -bool pmReadoutReadChunkMask(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, + return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_MASK, config); +} + +bool pmReadoutReadChunkMask(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap, pmConfig *config) { @@ -1100,5 +1104,5 @@ PS_ASSERT_INT_NONNEGATIVE(numScans, false); - return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_MASK, config); + return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_MASK, config); } @@ -1108,5 +1112,5 @@ PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_MASK, config); + return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_MASK, config); } @@ -1139,13 +1143,13 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config) +bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, pmConfig *config) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config); -} - -bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap, + return readoutMore(readout, fits, z, zMax, numScans, FPA_READ_TYPE_VARIANCE, config); +} + +bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int *zMax, int numScans, int overlap, pmConfig *config) { @@ -1155,5 +1159,5 @@ PS_ASSERT_INT_NONNEGATIVE(numScans, false); - return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config); + return readoutReadChunk(readout, fits, z, zMax, numScans, overlap, FPA_READ_TYPE_VARIANCE, config); } @@ -1163,5 +1167,5 @@ PS_ASSERT_FITS_NON_NULL(fits, false); - return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config); + return readoutReadChunk(readout, fits, z, NULL, 0, 0, FPA_READ_TYPE_VARIANCE, config); } Index: branches/pap/psModules/src/camera/pmFPARead.h =================================================================== --- branches/pap/psModules/src/camera/pmFPARead.h (revision 23948) +++ branches/pap/psModules/src/camera/pmFPARead.h (revision 25027) @@ -21,4 +21,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read pmConfig *config ///< Configuration @@ -31,4 +32,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read int overlap, ///< Overlap between consecutive reads @@ -98,4 +100,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read pmConfig *config ///< Configuration @@ -108,4 +111,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read int overlap, ///< Overlap between consecutive reads @@ -150,4 +154,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read pmConfig *config ///< Configuration @@ -160,4 +165,5 @@ psFits *fits, ///< FITS file from which to read int z, ///< Readout number/plane; zero-offset indexing + int *zMax, ///< Max plane number in this cell int numScans, ///< Number of scans (rows/cols) to read int overlap, ///< Overlap between consecutive reads Index: branches/pap/psModules/src/camera/pmFPAWrite.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAWrite.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPAWrite.c (revision 25027) @@ -160,9 +160,26 @@ psArray **imageArray = appropriateImageArray(hdu, type); // Array of images in the HDU + // XXX detect missing variance & mask images... + // Generate the HDU if needed --- this is required after a pmFPACopy, or similar, which does not // generate the HDU, but only copies the structure. - if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error."); - return false; + if (!blank && !hdu->blankPHU && !*imageArray) { + if (!pmHDUGenerateForCell(cell)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error."); + return false; + } + if (!*imageArray) { + if (type == FPA_WRITE_TYPE_IMAGE) { + psError(PS_ERR_UNKNOWN, false, "Expected to write an image, but it is missing...programming error?."); + return false; + } + if (type == FPA_WRITE_TYPE_MASK) { + psWarning("No mask image for this cell; skipping"); + } + if (type == FPA_WRITE_TYPE_VARIANCE) { + psWarning("No variance image for this cell; skipping"); + } + return true; + } } @@ -724,5 +741,5 @@ } if (!psFitsCompressionApply(fits, compress)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set FITS compression to NONE"); + psError(PS_ERR_UNKNOWN, false, "Unable to restore FITS compression"); psFree(compress); return false; Index: branches/pap/psModules/src/camera/pmFPAfileFringeIO.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileFringeIO.c (revision 25027) +++ branches/pap/psModules/src/camera/pmFPAfileFringeIO.c (revision 25027) @@ -0,0 +1,238 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmConfig.h" +//#include "pmConfigMask.h" +//#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +//#include "pmFPARead.h" +//#include "pmFPAWrite.h" +//#include "pmFPAMaskWeight.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAfileFringeIO.h" +#include "pmFPAfileFitsIO.h" +//#include "pmFPACopy.h" +//#include "pmFPAConstruct.h" +//#include "pmConceptsWrite.h" +#include "pmFringeStats.h" + +# define FRINGE_TABLE "FRINGE.MEASUREMENTS" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Reading Fringe data +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// given an already-opened fits file, read the table corresponding to the specified view +bool pmFPAviewReadFringes(const pmFPAview *view, pmFPAfile *file) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + pmFPA *fpa = file->fpa; // FPA of interest + psFits *fits = file->fits; // FITS file + + if (view->chip == -1) { + return pmFPAReadFringes(fpa, fits) > 0; + } + + if (view->cell == -1) { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + return pmChipReadFringes(chip, fits) > 0; + } + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + return pmCellReadFringes(cell, fits) > 0; +} + +int pmFPAReadFringes(pmFPA *fpa, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(fpa, 0); + PS_ASSERT_FITS_NON_NULL(fits, 0); + + int numRead = 0; // Number of reads + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + numRead += pmChipReadFringes(chip, fits); + } + + return numRead; +} + +int pmChipReadFringes(pmChip *chip, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(chip, 0); + PS_ASSERT_FITS_NON_NULL(fits, 0); + + int numRead = 0; // Number of reads + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + numRead += pmCellReadFringes(cell, fits); + } + + return numRead; +} + +int pmCellReadFringes(pmCell *cell, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(cell, 0); + PS_ASSERT_FITS_NON_NULL(fits, 0); + + const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + // Fringe Table Extension name + psString extname = NULL; + psStringAppend(&extname, "FRINGE_%s_%s", chipName, cellName); + + if (!psFitsMoveExtName(fits, extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname); + psFree(extname); + return 0; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname); + psFree(extname); + psFree(header); + return 0; + } + + psArray *table = psFitsReadTable(fits); // The table + if (!table) { + psError(PS_ERR_UNKNOWN, false, "Unable to add read the fringe table data from the file"); + psFree (header); + psFree (extname); + return 0; + } + + psArray *fringes = pmFringesParseTable(table, header); + psFree(table); + psFree(header); + + if (!fringes) { + psError(PS_ERR_UNKNOWN, false, "Unable to add parse the fringe table data"); + psFree(extname); + return 0; + } + + if (!psMetadataAdd(cell->analysis, PS_LIST_TAIL, FRINGE_TABLE, PS_DATA_ARRAY, "Fringes", fringes)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add fringe from extension %s to analysis metadata " + "for chip %s, cell %s\n", extname, chipName, cellName); + psFree(fringes); + psFree(extname); + return 0; + } + + psFree(fringes); // drop local reference + psFree(extname); + + return 1; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Writing fringe data +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// given an already-opened fits file, write the table corresponding to the specified view +bool pmFPAviewWriteFringes(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // FPA of interest + psFits *fits = file->fits; // FITS file + + if (view->chip == -1) { + return pmFPAWriteFringes(fits, fpa) > 0; + } + + if (view->cell == -1) { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + return pmChipWriteFringes(fits, chip) > 0; + } + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + return pmCellWriteFringes(fits, cell) > 0; +} + +int pmFPAWriteFringes(psFits *fits, const pmFPA *fpa) +{ + PS_ASSERT_PTR_NON_NULL(fpa, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + + int numWrite = 0; // Number of reads + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + numWrite += pmChipWriteFringes(fits, chip); + } + + return numWrite; +} + +int pmChipWriteFringes(psFits *fits, const pmChip *chip) +{ + PS_ASSERT_PTR_NON_NULL(chip, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + + int numWrite = 0; // Number of reads + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + numWrite += pmCellWriteFringes(fits, cell); + } + + return numWrite; +} + +int pmCellWriteFringes(psFits *fits, const pmCell *cell) +{ + PS_ASSERT_PTR_NON_NULL(cell, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + + const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + psArray *fringes = psMetadataLookupPtr(NULL, cell->analysis, FRINGE_TABLE); // The FITS table + if (!fringes) { + // We wrote everything we could find + return 0; + } + + psMetadata *header = psMetadataAlloc(); + psArray *table = pmFringesFormatTable(header, fringes); + if (!table) { + // We wrote everything we could find + psFree(header); + return 0; + } + + psString extname = NULL; // Extension name + psStringAppend(&extname, "FRINGE_%s_%s", chipName, cellName); + + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write table from chip %s, cell %s to extension %s\n", + chipName, cellName, extname); + psFree(extname); + psFree(header); + return 0; + } + + psFree(extname); + psFree(header); + return 1; +} + + Index: branches/pap/psModules/src/camera/pmFPAfileFringeIO.h =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileFringeIO.h (revision 25027) +++ branches/pap/psModules/src/camera/pmFPAfileFringeIO.h (revision 25027) @@ -0,0 +1,30 @@ +/* @file pmFPAfileFringeIO.h + * @brief Read & Write Fringe tables + * + * @author EAM, IfA + * @author PAP, IfA + * + * @version $Revision: 1.16 $ + * @date $Date: 2009-02-06 02:31:24 $ + * Copyright 2004-2005 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_FPA_FILE_FRINGE_IO_H +#define PM_FPA_FILE_FRINGE_IO_H + +/// @addtogroup Camera Camera Layout +/// @{ + +/// Read an fringes into the current view +bool pmFPAviewReadFringes(const pmFPAview *view, pmFPAfile *file); +int pmFPAReadFringes(pmFPA *fpa, psFits *fits); +int pmChipReadFringes(pmChip *chip, psFits *fits); +int pmCellReadFringes(pmCell *cell, psFits *fits); +bool pmFPAviewWriteFringes(const pmFPAview *view, pmFPAfile *file, pmConfig *config); +int pmFPAWriteFringes(psFits *fits, const pmFPA *fpa); +int pmChipWriteFringes(psFits *fits, const pmChip *chip); +int pmCellWriteFringes(psFits *fits, const pmCell *cell); + +/// @} + +# endif Index: branches/pap/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 23948) +++ branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 25027) @@ -23,4 +23,5 @@ #include "pmFPAWrite.h" #include "pmFPAfileFitsIO.h" +#include "pmFPAfileFringeIO.h" #include "pmSpan.h" #include "pmFootprint.h" @@ -192,5 +193,5 @@ status = pmFPAviewReadFitsImage(view, file, config); if (status) { - if (!pmFPAviewReadFitsTable(view, file, "FRINGE")) { + if (!pmFPAviewReadFringes(view, file)) { psError(PS_ERR_UNKNOWN, false, "Unable to read fringe data from %s.\n", file->filename); return false; @@ -451,5 +452,5 @@ status = pmFPAviewWriteFitsImage (view, file, config); if (status) { - if (!pmFPAviewWriteFitsTable(view, file, "FRINGE", config)) { + if (!pmFPAviewWriteFringes(view, file, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to write fringe data from %s.\n", file->filename); return false; Index: branches/pap/psModules/src/camera/pmHDU.c =================================================================== --- branches/pap/psModules/src/camera/pmHDU.c (revision 23948) +++ branches/pap/psModules/src/camera/pmHDU.c (revision 25027) @@ -93,12 +93,31 @@ } + psTrace("psModules.camera", 5, "Reading the header...\n"); + + // The header may already exist (e.g., from doing concept writing at the PHU level) so we need to be + // careful. We read into a separate container and copy that over the top of anything that's already read. + psMetadata *header = psFitsReadHeader(NULL, fits); + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); + return false; + } + if (!hdu->header) { - psTrace("psModules.camera", 5, "Reading the header...\n"); - hdu->header = psFitsReadHeader(hdu->header, fits); - if (! hdu->header) { - psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); - return false; - } - } + hdu->header = header; + return true; + } + + psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + const char *name = item->name; // Name of item + if (psMetadataLookup(hdu->header, name)) { + // It exists; clobber + psMetadataRemoveKey(hdu->header, name); + } + psMetadataAddItem(hdu->header, item, PS_LIST_TAIL, 0); + } + psFree(iter); + psFree(header); return true; Index: branches/pap/psModules/src/camera/pmHDUGenerate.c =================================================================== --- branches/pap/psModules/src/camera/pmHDUGenerate.c (revision 23948) +++ branches/pap/psModules/src/camera/pmHDUGenerate.c (revision 25027) @@ -390,7 +390,7 @@ if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) { // Nothing from which to create an HDU - psFree(cells); - psError(PS_ERR_IO, true, "Nothing from which to create an HDU\n"); - return false; + // psError(PS_ERR_IO, true, "Nothing from which to create an HDU\n"); + psWarning("Nothing from which to create an HDU, must be empty\n"); + return true; } @@ -504,6 +504,4 @@ psFree(iter); } - psFree(cells); - return true; } @@ -615,6 +613,7 @@ return true; } - - return generateHDU(hdu, cells); + bool status = generateHDU(hdu, cells); + psFree (cells); + return status; } case PM_FPA_LEVEL_CHIP: @@ -664,5 +663,7 @@ } - return generateHDU(hdu, cells); + bool status = generateHDU(hdu, cells); + psFree (cells); + return status; } case PM_FPA_LEVEL_FPA: Index: branches/pap/psModules/src/camera/pmReadoutFake.c =================================================================== --- branches/pap/psModules/src/camera/pmReadoutFake.c (revision 23948) +++ branches/pap/psModules/src/camera/pmReadoutFake.c (revision 25027) @@ -29,4 +29,6 @@ #define MAX_AXIS_RATIO 20.0 // Maximum axis ratio for PSF model #define SOURCE_MASK (PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources +#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ + PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models @@ -80,21 +82,4 @@ int numSources = sources->n; // Number of stars - - float flux0 = NAN; // Flux for central intensity of 1.0 - if (normalisePeak) { - pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0, - 1.0); // Fake model, with central intensity of 1.0 - psAssert(fakeModel, "failed to generate model: should this be an error or not?"); - - if (circularise && !circulariseModel(fakeModel)) { - psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); - psFree(fakeModel); - return false; - } - - flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0 - psFree(fakeModel); - } - for (int i = 0; i < numSources; i++) { pmSource *source = sources->data[i]; // Source of interest @@ -118,10 +103,25 @@ float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source + if (normalisePeak) { - flux /= flux0; + // Normalise flux + pmModel *normModel = pmModelFromPSFforXY(psf, x, y, 1.0); // Model for normalisation + if (!normModel || (normModel->flags & MODEL_MASK)) { + psFree(normModel); + continue; + } + if (circularise && !circulariseModel(normModel)) { + psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); + psFree(normModel); + return false; + } + + flux /= normModel->modelFlux(normModel->params); + psFree(normModel); } pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, flux); - if (!fakeModel) { + if (!fakeModel || (fakeModel->flags & MODEL_MASK)) { + psFree(fakeModel); continue; } Index: branches/pap/psModules/src/camera/pmReadoutStack.c =================================================================== --- branches/pap/psModules/src/camera/pmReadoutStack.c (revision 23948) +++ branches/pap/psModules/src/camera/pmReadoutStack.c (revision 25027) @@ -252,4 +252,7 @@ continue; } + if (!readout->process) { + continue; + } if (!readout->image) { psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %ld has NULL image.\n", i); @@ -260,7 +263,13 @@ pmCell *cell = readout->parent; // The parent cell bool mdok = true; // Status of MD lookup - psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section + psRegion *trimsec = cell ? psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC") : NULL; // Trim section if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) { - psWarning("CELL.TRIMSEC is not set for readout %ld --- ignored.\n", i); + psWarning("CELL.TRIMSEC is not set for readout %ld --- attempting to use image size.\n", i); + xSize = PS_MAX(xSize, readout->image->numCols); + ySize = PS_MAX(ySize, readout->image->numRows); + xMin = PS_MIN(xMin, 0); + xMax = PS_MAX(xMax, readout->image->numCols - 1); + yMin = PS_MIN(yMin, 0); + yMax = PS_MAX(yMax, readout->image->numRows - 1); } else { xSize = PS_MAX(xSize, trimsec->x1 - trimsec->x0); Index: branches/pap/psModules/src/concepts/pmConcepts.c =================================================================== --- branches/pap/psModules/src/concepts/pmConcepts.c (revision 23948) +++ branches/pap/psModules/src/concepts/pmConcepts.c (revision 25027) @@ -485,4 +485,28 @@ concept[length - 1] = '\0'; + // special variants: + if (!strcmp(concept, "FPA.DATE")) { + psTime *fpaTime = psMetadataLookupPtr(NULL, fpa->concepts, "FPA.TIME"); + if (!fpaTime) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Missing concept FPA.TIME needed for FPA.DATE"); + psFree(string); + return NULL; + } + psString dateTimeString = psTimeToISO(fpaTime); // String representation + psList *dateTime = psStringSplit(dateTimeString, "T", true); + psFree(dateTimeString); + psString dateString = psMemIncrRefCounter(psListGet(dateTime, PS_LIST_HEAD)); // The date string + psFree (dateTime); + + if (!psStringSubstitute(&string, dateString, "{FPA.DATE}")) { + psError(PS_ERR_UNKNOWN, false, "Unable to replace concept %s", concept); + psFree(string); + psFree(dateString); + return NULL; + } + psFree (dateString); + continue; + } + psTrace("psModules.concepts", 7, "Interpolating concept %s", concept); Index: branches/pap/psModules/src/concepts/pmConceptsRead.c =================================================================== --- branches/pap/psModules/src/concepts/pmConceptsRead.c (revision 23948) +++ branches/pap/psModules/src/concepts/pmConceptsRead.c (revision 25027) @@ -425,4 +425,5 @@ if (rows->n == 0) { psWarning("No rows returned from database query for concept %s --- ignored.", name); + psFree(rows); return NULL; } Index: branches/pap/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- branches/pap/psModules/src/concepts/pmConceptsStandard.c (revision 23948) +++ branches/pap/psModules/src/concepts/pmConceptsStandard.c (revision 25027) @@ -36,9 +36,9 @@ // Format type for time typedef enum { - TIME_FORMAT_YYYYMMDD, // Date stored in YYYY-MM-DD order (ISO-standard) - TIME_FORMAT_DDMMYYYY, // Date stored in DD-MM-YYYY order - TIME_FORMAT_MMDDYYYY, // Date stored in MM-DD-YYYY order - TIME_FORMAT_JD, // Date stored as JD - TIME_FORMAT_MJD, // Date stored as MJD + TIME_FORMAT_YYYYMMDD, // Date stored in YYYY-MM-DD order (ISO-standard) + TIME_FORMAT_DDMMYYYY, // Date stored in DD-MM-YYYY order + TIME_FORMAT_MMDDYYYY, // Date stored in MM-DD-YYYY order + TIME_FORMAT_JD, // Date stored as JD + TIME_FORMAT_MJD, // Date stored as MJD } conceptTimeFormatType; @@ -391,4 +391,5 @@ } + // FPA.RA and FPA.DEC psMetadataItem *p_pmConceptFormat_FPA_Coords(const psMetadataItem *concept, @@ -410,9 +411,8 @@ // How to interpret the coordinates bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMetadata(&mdok, - cameraFormat, - "FORMATS"); + psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS"); + bool sexagesimal = false; // Write sexagesimal format? if (mdok && formats) { - psString format = psMetadataLookupStr(&mdok,formats, concept->name); + psString format = psMetadataLookupStr(&mdok, formats, concept->name); if (mdok && strlen(format) > 0) { if (strcasecmp(format, "HOURS") == 0) { @@ -428,26 +428,36 @@ coords /= defaultCoordScaling(concept); } + + psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA + psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec + if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) { + sexagesimal = true; + } } else { coords /= defaultCoordScaling(concept); } - // We choose to write sexagesimal format - int big, medium; // Degrees and minutes - float small; // Seconds - bool negative = (coords < 0); // Are we working below zero? - coords = fabs(coords); - big = (int)abs(coords); - coords -= big; - medium = 60.0 * coords; - coords -= medium / 60.0; - small = 3600.0 * coords; - small = (float)((int)(small * 1000.0)) / 1000.0; - if (negative) { - big *= -1; - } - psString coordString = NULL; // String with the coordinates in sexagesimal format - psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small); - psMetadataItem *coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString); - psFree(coordString); + psMetadataItem *coordItem = NULL; // Item with coordinates, to return + if (sexagesimal) { + int big, medium; // Degrees and minutes + float small; // Seconds + bool negative = (coords < 0); // Are we working below zero? + coords = fabs(coords); + big = (int)abs(coords); + coords -= big; + medium = 60.0 * coords; + coords -= medium / 60.0; + small = 3600.0 * coords; + small = (float)((int)(small * 1000.0)) / 1000.0; + if (negative) { + big *= -1; + } + psString coordString = NULL; // String with the coordinates in sexagesimal format + psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small); + coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString); + psFree(coordString); + } else { + coordItem = psMetadataItemAllocF64(concept->name, concept->comment, coords); + } return coordItem; Index: branches/pap/psModules/src/config/pmConfig.c =================================================================== --- branches/pap/psModules/src/config/pmConfig.c (revision 23948) +++ branches/pap/psModules/src/config/pmConfig.c (revision 25027) @@ -243,4 +243,6 @@ unsigned int numBadLines = 0; struct stat filestat; + + pmErrorRegister(); psTrace("psModules.config", 3, "Loading %s configuration from file %s\n", @@ -1496,4 +1498,5 @@ } } + globfree(&globList); } psFree (words); Index: branches/pap/psModules/src/config/pmErrorCodes.dat =================================================================== --- branches/pap/psModules/src/config/pmErrorCodes.dat (revision 23948) +++ branches/pap/psModules/src/config/pmErrorCodes.dat (revision 25027) @@ -12,4 +12,5 @@ SKY Problem in sky STAMPS Unable to select stamps for PSF-matching +SMALL_AREA Small area for PSF-matching # these errors correspond to standard exit conditions ARGUMENTS Incorrect arguments Index: branches/pap/psModules/src/detrend/Makefile.am =================================================================== --- branches/pap/psModules/src/detrend/Makefile.am (revision 23948) +++ branches/pap/psModules/src/detrend/Makefile.am (revision 25027) @@ -16,5 +16,6 @@ pmShifts.c \ pmDark.c \ - pmRemnance.c + pmRemnance.c \ + pmPattern.c # pmSkySubtract.c @@ -33,5 +34,6 @@ pmShifts.h \ pmDark.h \ - pmRemnance.h + pmRemnance.h \ + pmPattern.h # pmSkySubtract.h Index: branches/pap/psModules/src/detrend/pmDark.c =================================================================== --- branches/pap/psModules/src/detrend/pmDark.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmDark.c (revision 25027) @@ -1,5 +1,10 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + #include #include #include +#include #include "psPolynomialMD.h" @@ -12,9 +17,10 @@ #include "pmReadoutStack.h" #include "pmDetrendThreads.h" - +#include "pmErrorCodes.h" #include "pmDark.h" #define PM_DARK_FITS_EXTNAME "PS_DARK" // FITS extension name for ordinates table #define PM_DARK_FITS_NAME "NAME" // Column name for concept name in ordinates table +#define PM_DARK_FITS_RULE "RULE" // Column name for concept rule in ordinates table #define PM_DARK_FITS_ORDER "ORDER" // Column name for polynomial order in ordinates table #define PM_DARK_FITS_SCALE "SCALE" // Column name for scaling option in ordinates table @@ -22,40 +28,24 @@ #define PM_DARK_FITS_MAX "MAX" // Column name for maximum value in ordinates table - - -// Look up the value of an ordinate in a readout -static bool ordinateLookup(float *value, // Value of ordinate, to return - bool *inRange, // is value within min : max range? - const char *name, // Name of ordinate (concept name) - bool scale, // Scale the value? - float min, float max, // Minimum and maximum values for scaling - const pmReadout *ro // Readout of interest - ) -{ - assert(value); - assert(name); - assert(ro); - - *inRange = true; - - pmCell *cell = ro->parent; // Parent cell - if (!cell) { - return false; - } +static bool ordinateParseConcept(double *value, const pmReadout *readout, const char *name) { + + *value = NAN; + + pmCell *cell = readout->parent; // Parent cell + psAssert(cell, "readout is missing cell \n"); + psMetadataItem *item = psMetadataLookup(cell->concepts, name); if (!item) { pmChip *chip = cell->parent; // Parent chip - if (!chip) { - return false; - } + psAssert(chip, "cell is missing chip \n"); + item = psMetadataLookup(chip->concepts, name); if (!item) { pmFPA *fpa = chip->parent; // Parent FPA - if (!fpa) { - return false; - } + psAssert(fpa, "chip is missing fpa \n"); + item = psMetadataLookup(fpa->concepts, name); if (!item) { - psWarning("Unable to find concept %s in readout", name); + psError(PM_ERR_CONFIG, true, "Unable to find concept %s in readout", name); return false; } @@ -63,9 +53,85 @@ } - *value = psMetadataItemParseF32(item); // Value of interest + *value = psMetadataItemParseF64(item); // Value of interest if (!isfinite(*value)) { psWarning("Non-finite value (%f) of concept %s in readout", *value, name); - return false; - } + } + return true; +} + +static bool ordinateParseRule(double *value, const pmReadout *readout, const char *name, const char *rule) { + + psAssert(name, "ordinate name is not defined"); + psAssert(rule, "ordinate rule is not defined"); + + *value = NAN; + + psArray *words = psStringSplitArray(rule, " ", false); + + // we should have a rule of the form (concept) OP (concept) OP (concept) ... + // for now, the only allowed OP is * (eventually, we can steal code from opihi for a better + // RPN parser). + + if (words->n % 2 == 0) { + psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + + for (int i = 1; i < words->n; i+=2) { + if (strcmp((char *)words->data[i], "*")) { + psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + } + + if (!ordinateParseConcept(value, readout, words->data[0])) { + psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + + double value2 = 0.0; + for (int i = 2; i < words->n; i+=2) { + if (!ordinateParseConcept(&value2, readout, words->data[i])) { + psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + *value *= value2; + } + psFree(words); + return true; +} + +// Look up the value of an ordinate in a readout +static bool ordinateLookup(double *value, // Value of ordinate, to return + bool *inRange, // is value within min : max range? + const char *name, // Name of ordinate (concept or abstract name) + const char *rule, // Rule for generating the value (if NULL use name as concept) + bool scale, // Scale the value? + float min, float max, // Minimum and maximum values for scaling + const pmReadout *readout // Readout of interest + ) +{ + assert(value); + assert(name); + assert(readout); + + *inRange = true; + + if (rule) { + if (!ordinateParseRule(value, readout, name, rule)) { + psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name); + return false; + } + } else { + if (!ordinateParseConcept(value, readout, name)) { + psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name); + return false; + } + } + if (scale) { if (*value < min || *value > max) { @@ -82,4 +148,5 @@ { psFree(ord->name); + psFree(ord->rule); return; } @@ -91,4 +158,5 @@ ord->name = psStringCopy(name); + ord->rule = NULL; ord->order = order; ord->scale = false; @@ -118,7 +186,11 @@ pmReadout *readout = inputs->data[i]; // Readout of interest - float normValue; // Normalisation value - if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { - psWarning("Unable to find value of %s for readout %d", normConcept, i); + double normValue; // Normalisation value + if (!ordinateLookup(&normValue, &inRange, normConcept, NULL, false, NAN, NAN, readout)) { + psError(PM_ERR_CONFIG, false, "problem finding concept %s for DARK.NORM", normConcept); + return false; + } + if (!isfinite(normValue)) { + psWarning("Unable to find acceptable value of %s for readout %d", normConcept, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; norm->data.F32[i] = NAN; @@ -140,5 +212,5 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information if (ord->order <= 0) { - psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order); + psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order); psFree(values); psFree(roMask); @@ -157,6 +229,11 @@ pmReadout *readout = inputs->data[j]; // Readout of interest - float value = NAN; // Value of ordinate - if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { + double value = NAN; // Value of ordinate + if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) { + psError(PM_ERR_CONFIG, false, "problem finding rule for DARK.ORDINATE %s", ord->name); + return false; + } + if (!isfinite(value)) { + psWarning("Unable to find acceptable value of DARK.ORDINATE %s for readout %d", ord->name, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; val->data.F32[i] = NAN; @@ -165,4 +242,5 @@ } if (!inRange) { + psWarning("Value of DARK.ORDINATE %s for readout %d is out of range", ord->name, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; val->data.F32[i] = NAN; @@ -309,4 +387,6 @@ return false; } + + pmDarkVisualInit(values); pmReadout *outReadout = output->readouts->data[0]; @@ -352,4 +432,8 @@ psVectorInit(poly->coeff, NAN); } + + pmDarkVisualPixelFit(pixels, mask); + pmDarkVisualPixelModel(poly, values); + for (int k = 0; k < poly->coeff->n; k++) { pmReadout *ro = output->readouts->data[k]; // Readout of interest @@ -447,16 +531,21 @@ for (int i = 0; i < numOrdinates; i++) { pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest - float value = NAN; // Value for ordinate - if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { - psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", ord->name); + double value = NAN; // Value for ordinate + if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) { + psError(PS_ERR_UNKNOWN, true, "Unable to find value for DARK.ORDINATE %s", ord->name); psFree(values); return false; } + if (!isfinite(value)) { + psError(PS_ERR_UNKNOWN, true, "Value for DARK.ORDINATE %s is NAN", ord->name); + psFree(values); + return false; + } values->data.F32[i] = value; } - float norm = NAN; // Normalisation value + double norm = NAN; // Normalisation value bool doNorm = false; // Do normalisation? if (normConcept && strlen(normConcept) > 0) { - if (!ordinateLookup(&norm, &inRange, normConcept, false, NAN, NAN, readout)) { + if (!ordinateLookup(&norm, &inRange, normConcept, NULL, false, NAN, NAN, readout)) { psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept); psFree(values); @@ -470,5 +559,5 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information if (ord->order <= 0) { - psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order); + psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order); psFree(values); psFree(orders); @@ -701,5 +790,13 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest psMetadata *row = psMetadataAlloc(); // FITS table row - psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name); + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "DARK.ORDINATE name", ord->name); + + // XXX write a dummy value if ord->rule == NULL? (eg, NONE) + if (ord->rule) { + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", ord->rule); + } else { + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", "NONE"); + } + psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order); psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale); @@ -878,4 +975,12 @@ ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX); + // load the ordinate rule; it is not an error if this field is missing or NULL + // a NULL rule means 'use the name as the single concept' + const char *rule = psMetadataLookupStr(&mdok, row, PM_DARK_FITS_RULE); + if (!rule || !strcasecmp(rule, "NONE")) { + ord->rule = NULL; + } else { + ord->rule = psStringCopy(rule); + } ordinates->data[i] = ord; } @@ -892,2 +997,152 @@ } +#if (HAVE_KAPA) +#include +#include "pmKapaPlots.h" +#include "pmVisual.h" + +static int nKapa = 0; +static int *kapa = NULL; +static bool plotFlag = true; +static psArray *xVectors = NULL; + +// this init function only gets the ordinates for the first readout... +bool pmDarkVisualInit(psArray *values) { + + if (!pmVisualIsVisual()) return true; + + // skip if we have already opened the windows (or if none are requested...) + if (nKapa) return true; + + // values has Ninput vectors with Norder elements; we need Norder vectors with Ninput elements... + int nOrders = 0; + for (int i = 0; i < values->n; i++) { + psVector *vect = values->data[i]; + if (!nOrders) { + nOrders = vect->n; + } else { + psAssert (nOrders == vect->n, "mismatch in order vector lengths"); + } + } + xVectors = psArrayAlloc(nOrders); + for (int i = 0; i < nOrders; i++) { + xVectors->data[i] = psVectorAlloc(values->n, PS_TYPE_F32); + } + + for (int i = 0; i < values->n; i++) { + psVector *vect = values->data[i]; + for (int j = 0; j < vect->n; j++) { + psVector *xVec = xVectors->data[j]; + xVec->data.F32[i] = vect->data.F32[j]; + } + } + + nKapa = nOrders; + + kapa = psAlloc(nKapa*sizeof(int)); + + for (int i = 0; i < nKapa; i++) { + kapa[i] = -1; + pmVisualInitWindow(&kapa[i], "ppmerge"); + } + return true; +} + +bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) { + + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + KapaInitGraph(&graphdata); + + psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n); + + psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32); + psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32); + + for (int i = 0; i < xVectors->n; i++) { + psVector *x = xVectors->data[i]; + + // generate vectors of the unmasked values + int nSub = 0; + for (int j = 0; j < pixels->n; j++) { + if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue; + xSub->data.F32[nSub] = x->data.F32[j]; + ySub->data.F32[nSub] = pixels->data.F32[j]; + nSub ++; + } + xSub->n = ySub->n = nSub; + + // plot the unmasked values + pmVisualScaleGraphdata (&graphdata, xSub, ySub, false); + KapaSetGraphData(kapa[i], &graphdata); + KapaSetLimits(kapa[i], &graphdata); + KapaClearPlots (kapa[i]); + + KapaSetFont (kapa[i], "courier", 14); + KapaBox (kapa[i], &graphdata); + KapaSendLabel (kapa[i], "ordinate", KAPA_LABEL_XM); + KapaSendLabel (kapa[i], "pixel values", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName("black"); + graphdata.style = 2; + graphdata.ptype = 2; + KapaPrepPlot (kapa[i], xSub->n, &graphdata); + KapaPlotVector(kapa[i], xSub->n, xSub->data.F32, "x"); + KapaPlotVector(kapa[i], xSub->n, ySub->data.F32, "y"); + } + pmVisualAskUser (&plotFlag); + return true; +} + +bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) { + + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + KapaInitGraph(&graphdata); + + psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n); + + psVector *yFit = psVectorAlloc(values->n, PS_TYPE_F32); + + for (int i = 0; i < values->n; i++) { + psVector *coord = values->data[i]; + yFit->data.F32[i] = psPolynomialMDEval (poly, coord); + } + + for (int i = 0; i < xVectors->n; i++) { + psVector *xFit = xVectors->data[i]; + + KapaGetGraphData(kapa[i], &graphdata); + graphdata.color = KapaColorByName("red"); + graphdata.style = 2; + graphdata.ptype = 7; + KapaPrepPlot (kapa[i], xFit->n, &graphdata); + KapaPlotVector(kapa[i], xFit->n, xFit->data.F32, "x"); + KapaPlotVector(kapa[i], xFit->n, yFit->data.F32, "y"); + } + pmVisualAskUser (&plotFlag); + return true; +} + +bool pmDarkVisualCleanup() { + + for (int i = 0; i < nKapa; i++) { + KapaClose(kapa[i]); + } + psFree (kapa); + psFree (xVectors); + return true; +} + +# else + +bool pmDarkVisualInit(psArray *values) { return true; } +bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) { return true; } +bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) { return true; } +bool pmDarkVisualCleanup() { return true; } + +# endif Index: branches/pap/psModules/src/detrend/pmDark.h =================================================================== --- branches/pap/psModules/src/detrend/pmDark.h (revision 23948) +++ branches/pap/psModules/src/detrend/pmDark.h (revision 25027) @@ -13,5 +13,6 @@ // An ordinate for fitting darks typedef struct { - psString name; // Name of concept to fit + psString name; // Name of ordinate + psString rule; // Rule for generating ordinate (math on concepts) int order; // Polynomial order to fit bool scale; // Rescale values? @@ -116,4 +117,8 @@ ); +bool pmDarkVisualInit(psArray *values); +bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask); +bool pmDarkVisualCleanup(); +bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values); #endif Index: branches/pap/psModules/src/detrend/pmDetrendDB.c =================================================================== --- branches/pap/psModules/src/detrend/pmDetrendDB.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmDetrendDB.c (revision 25027) @@ -105,4 +105,5 @@ DETREND_STRING_CASE(FRINGE); DETREND_STRING_CASE(ASTROM); + DETREND_STRING_CASE(NOISEMAP); default: return NULL; Index: branches/pap/psModules/src/detrend/pmDetrendDB.h =================================================================== --- branches/pap/psModules/src/detrend/pmDetrendDB.h (revision 23948) +++ branches/pap/psModules/src/detrend/pmDetrendDB.h (revision 25027) @@ -35,4 +35,5 @@ PM_DETREND_TYPE_BACKGROUND, PM_DETREND_TYPE_ASTROM, + PM_DETREND_TYPE_NOISEMAP, } pmDetrendType; Index: branches/pap/psModules/src/detrend/pmDetrendThreads.c =================================================================== --- branches/pap/psModules/src/detrend/pmDetrendThreads.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmDetrendThreads.c (revision 25027) @@ -60,4 +60,6 @@ } + // NOISEMAP : for now, not applied in the threaded loop + return true; } Index: branches/pap/psModules/src/detrend/pmFringeStats.c =================================================================== --- branches/pap/psModules/src/detrend/pmFringeStats.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmFringeStats.c (revision 25027) @@ -464,4 +464,8 @@ for (long i = 0; i < fringes->n; i++) { pmFringeStats *fringe = fringes->data[i]; // The fringe of interest + if (!fringe) { + psWarning ("skipping empty fringe stats for entry %ld -- video cell?\n", i); + continue; + } pmFringeRegions *regions = fringe->regions; // The fringe regions if (numPoints == 0) { @@ -496,4 +500,8 @@ for (long i = 0; i < fringes->n; i++) { pmFringeStats *fringe = fringes->data[i]; // The fringe of interest + if (!fringe) { + psWarning ("skipping empty fringe stats for entry %ld -- video cell?\n", i); + continue; + } pmFringeRegions *regions = fringe->regions; // The fringe regions // Copy the data over @@ -520,7 +528,7 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool pmFringesFormat(pmCell *cell, psMetadata *header, const psArray *fringes) -{ - PS_ASSERT_PTR_NON_NULL(cell, false); +psArray *pmFringesFormatTable(psMetadata *header, const psArray *fringes) +{ + PS_ASSERT_PTR_NON_NULL(header, false); PS_ASSERT_ARRAY_NON_NULL(fringes, false); @@ -531,5 +539,5 @@ if (stats->regions != regions) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Regions for fringe statistics are not identical.\n"); - return false; + return NULL; } } @@ -557,20 +565,8 @@ // Vectors: x, y, mask, f, df - psMetadata *scalars = psMemIncrRefCounter(header); // Metadata to hold the scalars; will be the header - if (!scalars) { - scalars = psMetadataAlloc(); - } - psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width", - regions->dX); - psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height", - regions->dY); - psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x", - regions->nX); - psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y", - regions->nY); - psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.HEADER", PS_DATA_METADATA, - "Header for fringe data", scalars); - psFree(scalars); - + psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width", regions->dX); + psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height", regions->dY); + psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x", regions->nX); + psMetadataAddS32(header, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y", regions->nY); psArray *table = psArrayAlloc(numRows); // The table @@ -605,27 +601,13 @@ } - psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE", PS_DATA_ARRAY, "Fringe data", table); - - return true; -} - - -psArray *pmFringesParse(pmCell *cell) -{ - PS_ASSERT_PTR_NON_NULL(cell, NULL); + return table; +} + +psArray *pmFringesParseTable(psArray *table, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(table, NULL); + PS_ASSERT_PTR_NON_NULL(header, NULL); bool mdok; // Status of MD lookup - psMetadata *header = psMetadataLookupMetadata(&mdok, cell->analysis, "FRINGE.HEADER"); // Header - if (!mdok || !header) { - psError(PS_ERR_UNKNOWN, false, "Unable to find header for fringe data.\n"); - return NULL; - } - - psArray *table = psMetadataLookupPtr(NULL, cell->analysis, "FRINGE"); // FITS table - if (!table) { - psError(PS_ERR_UNKNOWN, false, "Unable to find table for fringe data.\n"); - return NULL; - } - // Read the scalars from the header @@ -722,4 +704,55 @@ } +bool pmFringesFormat(pmCell *cell, psMetadata *inHeader, const psArray *fringes) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_ARRAY_NON_NULL(fringes, false); + + psMetadata *header = psMemIncrRefCounter(inHeader); // Metadata to hold the scalars; will be the header + if (!inHeader) { + inHeader = psMetadataAlloc(); + } + + psArray *table = pmFringesFormatTable(header, fringes); + if (!table) { + psError (PS_ERR_UNKNOWN, false, "unable to generate table from fringes"); + psFree(header); + return NULL; + } + + psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.HEADER", PS_DATA_METADATA, "Header for fringe data", header); + psFree(header); + + psMetadataAdd(cell->analysis, PS_LIST_TAIL, "FRINGE.TABLE", PS_DATA_ARRAY, "Fringe data", table); + psFree(table); + + return true; +} + +psArray *pmFringesParse(pmCell *cell) +{ + PS_ASSERT_PTR_NON_NULL(cell, NULL); + + bool mdok; // Status of MD lookup + psMetadata *header = psMetadataLookupMetadata(&mdok, cell->analysis, "FRINGE.HEADER"); // Header + if (!mdok || !header) { + psError(PS_ERR_UNKNOWN, false, "Unable to find header for fringe data.\n"); + return NULL; + } + + psArray *table = psMetadataLookupPtr(NULL, cell->analysis, "FRINGE.TABLE"); // FITS table + if (!table) { + psError(PS_ERR_UNKNOWN, false, "Unable to find table for fringe data.\n"); + return NULL; + } + + psArray *fringes = pmFringesParseTable(table, header); + if (!fringes) { + psError(PS_ERR_UNKNOWN, false, "Unable to add parse the fringe table data"); + return NULL; + } + + return fringes; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -821,6 +854,8 @@ // Solve the least-squares equation - if (! psMatrixGJSolve(A, B)) { - psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); + if (!psMatrixGJSolve(A, B)) { + psLogMsg("psModules.detrend", PS_LOG_INFO, "Could not solve linear equations. Returning NULL.\n"); + psFree(A); + psFree(B); return false; } @@ -882,5 +917,8 @@ psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics - psVectorStats(stats, diffs, NULL, mask, 1); + if (!psVectorStats(stats, diffs, NULL, mask, 1)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return 0; + } float middle = stats->sampleMedian; // The middle of the distribution float thresh = rej * 0.74 * (stats->sampleUQ - stats->sampleLQ); // The rejection threshold @@ -969,9 +1007,15 @@ // Get rid of the extreme outliers by assuming most of the points are somewhat clustered - psVectorStats(median, science->f, NULL, NULL, 0); + if (!psVectorStats(median, science->f, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return NULL; + } scale->coeff->data.F32[0] = median->sampleMedian; for (int i = 0; i < fringes->n; i++) { pmFringeStats *fringe = fringes->data[i]; // The fringe of interest - psVectorStats(median, fringe->f, NULL, NULL, 0); + if (!psVectorStats(median, fringe->f, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return NULL; + } scale->coeff->data.F32[0] -= median->sampleMedian; scale->coeff->data.F32[i] = 0.0; Index: branches/pap/psModules/src/detrend/pmFringeStats.h =================================================================== --- branches/pap/psModules/src/detrend/pmFringeStats.h (revision 23948) +++ branches/pap/psModules/src/detrend/pmFringeStats.h (revision 25027) @@ -142,11 +142,27 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -/// Write an array of fringes measurements to a FITS table. -/// -/// Writes an array of fringe measurements for a cell as a FITS table in the analysis metadata. The array of -/// fringe statistics must all use the same fringe regions (or there is no point in storing them all -/// together). The header is supplemented with scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY, -/// PSFRNGNX, PSFRNGNY) from the fringe regions, while the fringe coordinates and mask are written as a FITS -/// table (as x, y, mask, f, df; f and df are vectors). +/// Convert an array of fringes measurements to a psArray suitable for writing as a FITS table. +/// +/// Converts an array of fringe measurements for a cell into the corresponding rows of a FITS +/// table (array of psMetadata). The array of fringe statistics must all use the same fringe +/// regions (or there is no point in storing them all together). The header is supplemented +/// with scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY, PSFRNGNX, PSFRNGNY) from the +/// fringe regions, while the fringe coordinates and mask are written as a FITS table rows (as +/// x, y, mask, f, df; f and df are vectors). Use psFitsTableWrite to save the resulting rows +/// to disk. +psArray *pmFringesFormatTable(psMetadata *header, const psArray *fringes); + + +/// Parses an array of fringes measurements from a FITS table. +/// +/// The fringes for the cell are read from the FITS table (array of psMetadata rows). The +/// table provides the region and the (possibly multiple) fringe statistics for that region. +/// The supplied header defines the scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY, +/// PSFRNGNX, PSFRNGNY) +psArray *pmFringesParseTable(psArray *table, psMetadata *header); + +/// Deprecated Function: converts the fringes to a FITS table (array of psMetadata) and saves +/// them on the cell->analysis; scalar values are dX, dY, nX and nY are written to the header +/// (as PSFRNGDX, PSFRNGDY, PSFRNGNX, PSFRNGNY) bool pmFringesFormat(pmCell *cell, ///< Cell for which to write psMetadata *header, ///< Header, or NULL @@ -154,12 +170,8 @@ ); -/// Parses an array of fringes measurements from a FITS table. -/// -/// The fringes for the cell are read from the FITS table in the analysis metadata. The table provides the -/// region and the (possibly multiple) fringe statistics for that region. The current extension is used if -/// the extension name is not provided. +/// Deprecated Function: pulls a the header and FITS table (array of psMetadata) representing +/// the fringes from the cell->analysis and converts to fringe measurements psArray *pmFringesParse(pmCell *cell ///< Cell for which to read fringes ); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// Index: branches/pap/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- branches/pap/psModules/src/detrend/pmMaskBadPixels.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmMaskBadPixels.c (revision 25027) @@ -81,6 +81,6 @@ -bool pmMaskFlagSuspectPixels(pmReadout *output, const pmReadout *readout, float median, float stdev, - float rej, psImageMaskType maskVal) +bool pmMaskFlagSuspectPixelsBySigma(pmReadout *output, const pmReadout *readout, float median, float stdev, + float rej, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -115,5 +115,6 @@ // If we get down here and the statistics are missing, then we should go and mask the entire image psWarning("Missing statistics --- flagging entire image as suspect."); - return (psImage*)psBinaryOp(suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp (suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_F32)); + return true; } @@ -128,4 +129,59 @@ for (int x = 0; x < image->numCols; x++) { if (fabs((image->data.F32[y][x] - median) / stdev) < rej) continue; + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) continue; + suspect->data.F32[y][x] += 1.0; + } + } + psFree(suspect); // Drop reference + + psMetadataItem *numItem = psMetadataLookup(output->analysis, PM_MASK_ANALYSIS_NUM); // Item with number + assert(numItem); + numItem->data.S32++; + + return true; +} + +bool pmMaskFlagSuspectPixelsByValue(pmReadout *output, const pmReadout *readout, + float min, float max, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_IMAGE_NON_NULL(readout->image, false); + PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false); + PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false); + if (readout->mask) { + PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false); + PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_IMAGE_MASK, false); + } + PS_ASSERT_PTR_NON_NULL(output, false); + + bool mdok; // Status of MD lookup + psImage *suspect = psMetadataLookupPtr(&mdok, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img + if (suspect) { + PS_ASSERT_IMAGE_NON_EMPTY(suspect, false); + PS_ASSERT_IMAGE_TYPE(suspect, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, suspect, false); + psMemIncrRefCounter(suspect); + } else { + suspect = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); + psImageInit(suspect, 0); + psMetadataAddImage(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_SUSPECT, PS_META_REPLACE, + "Suspect pixels", suspect); + psMetadataAddS32(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_NUM, PS_META_REPLACE, + "Number of input images", 0); + } + + psImage *image = readout->image; // Image of interest + psImage *mask = readout->mask; // Corresponding mask + + psTrace ("psModules.detrend", 3, "suspect: < %f or > %f\n", min, max); + + // XXX this loop could result in pixels with suspect = 0.0 but no valid input pixels (all + // masked). need to track the number of good as well as suspect pixels? + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + bool above = image->data.F32[y][x] > max; + bool below = image->data.F32[y][x] < min; + if (!above && !below) continue; if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) continue; suspect->data.F32[y][x] += 1.0; Index: branches/pap/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- branches/pap/psModules/src/detrend/pmMaskBadPixels.h (revision 23948) +++ branches/pap/psModules/src/detrend/pmMaskBadPixels.h (revision 25027) @@ -50,10 +50,24 @@ /// high value in the suspect pixels image, allowing them to be identified. The suspect pixels /// image is of type S32. The relevant median and standard deviation must be supplied in the -/// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEVe -bool pmMaskFlagSuspectPixels(pmReadout *output, ///< Output readout, optionally with suspect pixels image +/// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEV +bool pmMaskFlagSuspectPixelsBySigma(pmReadout *output, ///< Output readout, optionally with suspect pixels image const pmReadout *readout, ///< Readout to inspect float median, ///< Image median float stdev, ///< Image standard deviation float rej, ///< Rejection threshold (standard deviations) + psImageMaskType maskVal ///< Mask value for statistics + ); + +/// Find out-of-range pixels and flag them +/// +/// Pixels great > max or < min have the corresponding pixel in the "suspect pixels" image +/// incremented. After accumulating over a suitable sample of images, bad pixels should have a +/// high value in the suspect pixels image, allowing them to be identified. The suspect pixels +/// image is of type S32. The relevant median and standard deviation must be supplied in the +/// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEV +bool pmMaskFlagSuspectPixelsByValue(pmReadout *output, ///< Output readout, optionally with suspect pixels image + const pmReadout *readout, ///< Readout to inspect + float min, ///< Image min acceptable value + float max, ///< Image max acceptable value psImageMaskType maskVal ///< Mask value for statistics ); Index: branches/pap/psModules/src/detrend/pmOverscan.c =================================================================== --- branches/pap/psModules/src/detrend/pmOverscan.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmOverscan.c (revision 25027) @@ -74,5 +74,8 @@ mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] - psVectorStats(myStats, values, NULL, NULL, 0); + if (!psVectorStats(myStats, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } reduced->data.F32[i] = psStatsGetValue(myStats, statistic); } else if (overscanOpts->fitType == PM_FIT_NONE) { @@ -275,5 +278,8 @@ psFree(iter); - (void)psVectorStats(stats, pixels, NULL, NULL, 0); + if (!psVectorStats(stats, pixels, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } psFree(pixels); double reduced = psStatsGetValue(stats, statistic); // Result of statistics @@ -349,5 +355,8 @@ psString comment = NULL; // Comment to add psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (vectorStats, reduced, NULL, NULL, 0); + if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); @@ -412,5 +421,8 @@ psString comment = NULL; // Comment to add psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (vectorStats, reduced, NULL, NULL, 0); + if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); Index: branches/pap/psModules/src/detrend/pmPattern.c =================================================================== --- branches/pap/psModules/src/detrend/pmPattern.c (revision 25027) +++ branches/pap/psModules/src/detrend/pmPattern.c (revision 25027) @@ -0,0 +1,127 @@ +#include +#include + +#include "pmFPA.h" + +#include "pmPattern.h" + + +// Mask a row as bad +static void patternMaskRow(pmReadout *ro, // Readout to mask + int y, // Row to mask + psImageMaskType bad // Mask value to give + ) +{ + psImage *image = ro->image; // Image to mask + psAssert(image, "Require image"); + psAssert(y < image->numRows, "Row not in image"); + + int numCols = image->numCols; // Size of image + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] = NAN; + } + if (ro->mask) { + psImage *mask = ro->mask; // Mask image to mask + for (int x = 0; x < numCols; x++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= bad; + } + } + return; +} + +bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PS_ASSERT_INT_NONNEGATIVE(order, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + psImage *image = ro->image; // Image to correct + psImage *mask = ro->mask; // Mask for image + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { + psWarning("Unable to calculate statistics on readout; masking entire readout."); + psErrorClear(); + psFree(stats); + psFree(rng); + psImageInit(image, NAN); + if (mask) { + psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); + } + if (ro->variance) { + psImageInit(image, NAN); + } + return true; + } + float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data + float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data + psFree(stats); + psFree(rng); + + // Indices are distributed [-1:1) + psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fitting + float norm = 2.0 / (float)numCols; // Normalisation for indices + for (int x = 0; x < numCols; x++) { + indices->data.F32[x] = x * norm - 1.0; + } + + psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics + clip->clipIter = iter; + clip->clipSigma = rej; + psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit + psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit + + for (int y = 0; y < numRows; y++) { + psVectorInit(clipMask, 0); + data = psImageRow(data, image, y); + int num = 0; // Number of good pixels + for (int x = 0; x < numCols; x++) { + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || + data->data.F32[x] < lower || data->data.F32[x] > upper) { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; + } else { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0; + num++; + } + } + if (num < order + 1) { + // Not enough points to fit + patternMaskRow(ro, y, maskBad); + continue; + } + if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) { + psWarning("Unable to fit polynomial to row %d", y); + psErrorClear(); + patternMaskRow(ro, y, maskBad); + continue; + } + psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector + if (!solution) { + psWarning("Unable to evaluate polynomial for row %d", y); + psErrorClear(); + patternMaskRow(ro, y, maskBad); + continue; + } + + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] -= solution->data.F32[x]; + } + psFree(solution); + } + + psFree(indices); + psFree(clip); + psFree(clipMask); + psFree(poly); + psFree(data); + + return true; +} Index: branches/pap/psModules/src/detrend/pmPattern.h =================================================================== --- branches/pap/psModules/src/detrend/pmPattern.h (revision 25027) +++ branches/pap/psModules/src/detrend/pmPattern.h (revision 25027) @@ -0,0 +1,37 @@ +/* @file pmPattern.h + * @brief Fit and remove pattern noise + * + * @author Paul Price, IfA + * + * @version $Revision: 1.16 $ + * @date $Date: 2009-02-12 19:25:52 $ + * Copyright 2004-2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_PATTERN_H +#define PM_PATTERN_H + +#include +#include + +/// @addtogroup detrend Detrend Creation and Application +/// @{ + +/// Fit and remove pattern noise over rows +bool pmPatternRow( + pmReadout *ro, ///< Readout to correct + int order, ///< Polynomial order + int iter, ///< Number of clipping iterations + float rej, ///< Rejection threshold for clipping + float thresh, ///< Threshold for ignoring pixels + psStatsOptions clipMean, ///< Statistic to use for mean + psStatsOptions clipStdev, ///< Statistic to use for standard deviation + psImageMaskType maskVal, ///< Mask value to use + psImageMaskType maskBad ///< Mask value to give bad pixels + ); + +/// @} +#endif + + + Index: branches/pap/psModules/src/detrend/pmRemnance.c =================================================================== --- branches/pap/psModules/src/detrend/pmRemnance.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmRemnance.c (revision 25027) @@ -66,6 +66,8 @@ values->n = numValues; if (!psVectorStats(stats, values, NULL, NULL, 0)) { - // Can't do anything about it - psErrorClear(); + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + if (isnan(stats->sampleMedian)) { maxMask = max; continue; Index: branches/pap/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 25027) @@ -350,10 +350,17 @@ psStats *rawStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (rawStats, counts, NULL, NULL, 0); - psVectorStats (resStats, resid, NULL, NULL, 0); + if (!psVectorStats (rawStats, counts, NULL, NULL, 0)) { + 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; + } // XXX temporary hard-wired minimum stdev improvement factor psTrace("psModules.detrend", 3, "raw scatter %f vs res scatter %f\n", rawStats->sampleStdev, resStats->sampleStdev); if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false; + if (isnan(rawStats->sampleStdev) || isnan(resStats->sampleStdev)) corr->valid = false; psFree (rawStats); Index: branches/pap/psModules/src/detrend/pmSkySubtract.c =================================================================== --- branches/pap/psModules/src/detrend/pmSkySubtract.c (revision 23948) +++ branches/pap/psModules/src/detrend/pmSkySubtract.c (revision 25027) @@ -410,9 +410,9 @@ // XXX: How do we know if these matrix operations were successful? // - Aout = psMatrixLUD(Aout, &outPerm, A); + Aout = psMatrixLUDecomposition(Aout, &outPerm, A); PS_ASSERT_IMAGE_NON_NULL(Aout, NULL); PS_ASSERT_IMAGE_NON_EMPTY(Aout, NULL); psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64); - psMatrixLUSolve(C, Aout, B, outPerm); + psMatrixLUSolution(C, Aout, B, outPerm); // Index: branches/pap/psModules/src/extras/Makefile.am =================================================================== --- branches/pap/psModules/src/extras/Makefile.am (revision 23948) +++ branches/pap/psModules/src/extras/Makefile.am (revision 25027) @@ -8,5 +8,6 @@ psIOBuffer.c \ pmKapaPlots.c \ - pmVisual.c + pmVisual.c \ + ippStages.c pkginclude_HEADERS = \ @@ -15,5 +16,6 @@ psIOBuffer.h \ pmKapaPlots.h \ - pmVisual.h + pmVisual.h \ + ippStages.h CLEANFILES = *~ Index: branches/pap/psModules/src/extras/ippStages.c =================================================================== --- branches/pap/psModules/src/extras/ippStages.c (revision 25027) +++ branches/pap/psModules/src/extras/ippStages.c (revision 25027) @@ -0,0 +1,31 @@ +/** Functions that describe the various stages of the IPP + * @author Bill Sweeney, IfA + **/ + +#include +#include "ippStages.h" +#include + +ippStage ippStringToStage(const psString stageString) +{ + PS_ASSERT_PTR_NON_NULL(stageString, IPP_STAGE_NONE); + + if (!strcmp(stageString, "raw")) { + return IPP_STAGE_RAW; + } else if (!strcmp(stageString, "chip")) { + return IPP_STAGE_CHIP; + } else if (!strcmp(stageString, "camera")) { + return IPP_STAGE_CAMERA; + } else if (!strcmp(stageString, "warp")) { + return IPP_STAGE_WARP; + } else if (!strcmp(stageString, "fake")) { + return IPP_STAGE_FAKE; + } else if (!strcmp(stageString, "diff")) { + return IPP_STAGE_DIFF; + } else if (!strcmp(stageString, "stack")) { + return IPP_STAGE_STACK; + } else { + psError(PS_ERR_PROGRAMMING, true, "%s is not a valid IPP stage", stageString); + return IPP_STAGE_NONE; + } +} Index: branches/pap/psModules/src/extras/ippStages.h =================================================================== --- branches/pap/psModules/src/extras/ippStages.h (revision 25027) +++ branches/pap/psModules/src/extras/ippStages.h (revision 25027) @@ -0,0 +1,27 @@ +/* @file psVisual.h + * @brief some macro defintions for the stages of the pipeline + * @author Bill Sweeney, IfA + * + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef IPP_STAGES_H +#define IPP_STAGES_H + +typedef enum { + IPP_STAGE_NONE = -1, + IPP_STAGE_RAW = 0, + IPP_STAGE_CHIP, + IPP_STAGE_CAMERA, + IPP_STAGE_FAKE, + IPP_STAGE_WARP, + IPP_STAGE_DIFF, + IPP_STAGE_STACK, +} ippStage; + +/** return the ippStage represented by a string + * @return the corresponding value or IPP_STAGE_NONE if invalid + */ +ippStage ippStringToStage(const psString stageString); + +#endif Index: branches/pap/psModules/src/extras/pmVisual.c =================================================================== --- branches/pap/psModules/src/extras/pmVisual.c (revision 23948) +++ branches/pap/psModules/src/extras/pmVisual.c (revision 25027) @@ -281,6 +281,12 @@ psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - psVectorStats (statsX, xVec, NULL, NULL, 0); - psVectorStats (statsY, yVec, NULL, NULL, 0); + if (!psVectorStats (statsX, xVec, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + if (!psVectorStats (statsY, yVec, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } float xhi = statsX->sampleMedian + 3 *statsX->sampleStdev; Index: branches/pap/psModules/src/imcombine/pmImageCombine.c =================================================================== --- branches/pap/psModules/src/imcombine/pmImageCombine.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmImageCombine.c (revision 25027) @@ -132,4 +132,8 @@ // Combine all the pixels, using the specified stat. if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + if (isnan(stats->sampleMean)) { combine->data.F32[y][x] = NAN; psFree(buffer); @@ -137,5 +141,5 @@ } float combinedPixel = stats->sampleMean; // Value of the combination - + if (iter == 0) { // Save the value produced with no rejection, since it may be useful later @@ -364,5 +368,7 @@ // Get the median psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - psVectorStats(stats, pixels, NULL, mask, 1); + if (!psVectorStats(stats, pixels, NULL, mask, 1)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + } float median = stats->sampleMedian; psFree(stats); Index: branches/pap/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/pap/psModules/src/imcombine/pmPSFEnvelope.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmPSFEnvelope.c (revision 25027) @@ -33,5 +33,5 @@ -//#define TESTING // Enable test output +// #define TESTING // Enable test output #define PEAK_FLUX 1.0e4 // Peak flux for each source #define SKY_VALUE 0.0e0 // Sky value for fake image @@ -40,4 +40,6 @@ #define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF #define SOURCE_FIT_ITERATIONS 100 // Number of iterations for source fitting +#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ + PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models @@ -112,5 +114,7 @@ psImageInit(envelope, SKY_VALUE); pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout - float maxRadius = 0.0; // Maximum radius of sources + float maxRadius = 0.0; // Maximum radius for sources + psVector *numbers = psVectorAlloc(numFakes, PS_TYPE_S32); // Number of detections for each source + psVectorInit(numbers, 0); for (int i = 0; i < inputs->n; i++) { pmPSF *psf = inputs->data[i]; // PSF of interest @@ -127,4 +131,5 @@ psFree(xOffset); psFree(fakes); + psFree(numbers); psf->residuals = resid; return NULL; @@ -140,4 +145,7 @@ double flux = fakeRO->image->data.F32[(int)y][(int)x]; + if (!isfinite(flux) || flux < 0) { + continue; + } float norm = PEAK_FLUX / flux; // Normalisation for source psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region @@ -151,10 +159,17 @@ // Get the radius pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source - psAssert (model, "failed to generate model: should this be an error or not?"); + if (!model || (model->flags & MODEL_MASK)) { + continue; + } float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source + if (srcRadius == 0) { + continue; + } if (srcRadius > maxRadius) { maxRadius = srcRadius; } - psFree(model); + + // If we got this far, the source is decent + numbers->data.S32[j]++; } @@ -175,8 +190,4 @@ psFree(fakeRO); - if (maxRadius > radius) { - maxRadius = radius; - } - #ifdef TESTING { @@ -190,4 +201,5 @@ // Put the fake sources onto a full-size image + psArray *goodFakes = psArrayAllocEmpty(numFakes); // Good fake sources pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -195,4 +207,8 @@ for (int i = 0; i < numFakes; i++) { pmSource *source = fakes->data[i]; // Fake source + if (numbers->data.S32[i] > 0) { + psArrayAdd(goodFakes, goodFakes->n, source); + } + // Position of source on fake image int xFake = source->peak->x + xOffset->data.S32[i]; @@ -213,4 +229,5 @@ psFree(yOffset); psFree(fakes); + psFree(numbers); return NULL; } @@ -220,4 +237,16 @@ psFree(yOffset); psFree(envelope); + psFree(numbers); + + psFree(fakes); + fakes = goodFakes; + numFakes = fakes->n; + + if (numFakes == 0.0) { + psError(PS_ERR_UNKNOWN, false, "No fake sources are suitable for PSF fitting."); + psFree(fakes); + psFree(readout); + return false; + } // XXX Setting the variance seems to be an art @@ -232,4 +261,8 @@ psImageInit(readout->mask, 0); + if (maxRadius > radius) { + maxRadius = radius; + } + #ifdef TESTING { @@ -243,4 +276,5 @@ // Reset the sources to point to the new pixels, and measure the moments in preparation for PSF fitting + int numMoments = 0; // Number of moments measured for (int i = 0; i < numFakes; i++) { pmSource *source = fakes->data[i]; // Fake source @@ -257,5 +291,5 @@ source->maskObj = NULL; - if (!pmSourceDefinePixels(source, readout, x, y, maxRadius)) { + if (!pmSourceDefinePixels(source, readout, x, y, radius)) { psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source."); psFree(readout); @@ -264,10 +298,18 @@ } - if (!pmSourceMoments(source, maxRadius)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source."); - psFree(readout); - psFree(fakes); - return NULL; - } + // measure the source moments: tophat windowing, no pixel S/N cutoff + if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) { + // Can't do anything about it; limp along as best we can + psErrorClear(); + continue; + } + numMoments++; + } + + if (numMoments == 0) { + psError(PS_ERR_UNKNOWN, true, "Unable to measure moments for sources."); + psFree(fakes); + psFree(readout); + return NULL; } @@ -286,4 +328,5 @@ options->psfFieldXo = 0; options->psfFieldYo = 0; + options->chiFluxTrend = false; // All sources have similar flux, so fitting a trend often fails pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); Index: branches/pap/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- branches/pap/psModules/src/imcombine/pmReadoutCombine.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmReadoutCombine.c (revision 25027) @@ -173,4 +173,6 @@ for (int i = 0; i < inputs->n; i++) { pmReadout *readout = inputs->data[i]; // Readout of interest + psAssert(readout, "readout was not defined"); + if (!readout->process) continue; if (!readout->image) { psError(PS_ERR_UNEXPECTED_NULL, true, "Image data is missing for image %d.\n", i); @@ -271,4 +273,5 @@ for (int r = 0; r < inputs->n; r++) { pmReadout *readout = inputs->data[r]; // Input readout + if (!readout->process) continue; psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r, minInputCols - readout->col0, maxInputCols - readout->col0, @@ -277,4 +280,7 @@ } #endif + + // set up windows for visualization (if selected) + pmReadoutCombineVisualInit(); // Dereference output products @@ -308,4 +314,8 @@ for (int r = 0; r < inputs->n; r++) { pmReadout *readout = inputs->data[r]; // Input readout + if (!readout->process) { + maskData[r] = 1; + continue; + } int yIn = i - readout->row0; // y position on input readout int xIn = j - readout->col0; // x position on input readout @@ -384,23 +394,31 @@ // Combination if (!psVectorStats(stats, pixels, errors, mask, 1)) { - // Can't do much about it, but it's not worth worrying about - psErrorClear(); + psError(PS_ERR_UNKNOWN, false, "error in pixel stats"); + return false; + } + + outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); + + if (!isfinite(outputImage[yOut][xOut])) { + pmReadoutCombineVisualPixels(pixels, mask, outputImage[yOut][xOut]); + } + + if (isnan(outputImage[yOut][xOut])) { outputImage[yOut][xOut] = NAN; outputMask[yOut][xOut] = params->blank; + sigma->data.F32[yOut][xOut] = NAN; if (params->variances) { outputVariance[yOut][xOut] = NAN; } - sigma->data.F32[yOut][xOut] = NAN; - } else { - outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); - outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank; - if (params->variances) { - float stdev = psStatsGetValue(stats, combineStdev); - outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance - // XXXX this is not the correct formal error. - // also, the weighted mean is not obviously the correct thing here - } - sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev); - } + continue; + } + outputMask[yOut][xOut] = 0; + sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev); + if (params->variances) { + float stdev = psStatsGetValue(stats, combineStdev); + outputVariance[yOut][xOut] = PS_SQR(stdev); // Variance + // XXXX this is not the correct formal error. + // also, the weighted mean is not obviously the correct thing here + } } } @@ -423,2 +441,106 @@ } +#if (HAVE_KAPA) +#include +#include "pmKapaPlots.h" +#include "pmVisual.h" + +static int kapa = -1; +static bool plotFlag = true; + +// this init function only gets the ordinates for the first readout... +bool pmReadoutCombineVisualInit(void) { + + if (!pmVisualIsVisual()) return true; + + // skip if we have already opened the windows (or if none are requested...) + if (kapa != -1) return true; + + pmVisualInitWindow(&kapa, "ppmerge"); + return true; +} + +bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) { + + Graphdata graphdata; + float xline[2], yline[2]; + + if (!pmVisualIsVisual()) return true; + + if (!plotFlag) return true; + + KapaInitGraph(&graphdata); + + psVector *xAll = psVectorAlloc(pixels->n, PS_TYPE_F32); + psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32); + psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32); + + // generate vectors of the unmasked values + int nSub = 0; + for (int j = 0; j < pixels->n; j++) { + xAll->data.F32[j] = j; + if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue; + xSub->data.F32[nSub] = j; + ySub->data.F32[nSub] = pixels->data.F32[j]; + nSub ++; + } + xSub->n = ySub->n = nSub; + xAll->n = pixels->n; + + xline[0] = 0; + xline[1] = pixels->n; + yline[0] = mean; + yline[1] = mean; + + // plot the unmasked values + pmVisualScaleGraphdata (&graphdata, xAll, pixels, false); + KapaSetGraphData(kapa, &graphdata); + KapaSetLimits(kapa, &graphdata); + KapaClearPlots (kapa); + + KapaSetFont (kapa, "courier", 14); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "ordinate", KAPA_LABEL_XM); + KapaSendLabel (kapa, "pixel values", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName("black"); + graphdata.style = 2; + graphdata.ptype = 2; + KapaPrepPlot (kapa, xSub->n, &graphdata); + KapaPlotVector(kapa, xSub->n, xSub->data.F32, "x"); + KapaPlotVector(kapa, xSub->n, ySub->data.F32, "y"); + + graphdata.color = KapaColorByName("red"); + graphdata.style = 2; + graphdata.ptype = 7; + KapaPrepPlot (kapa, xAll->n, &graphdata); + KapaPlotVector(kapa, xAll->n, xAll->data.F32, "x"); + KapaPlotVector(kapa, xAll->n, pixels->data.F32, "y"); + + graphdata.color = KapaColorByName("blue"); + graphdata.style = 0; + graphdata.ptype = 7; + KapaPrepPlot (kapa, 2, &graphdata); + KapaPlotVector(kapa, 2, xline, "x"); + KapaPlotVector(kapa, 2, yline, "y"); + + pmVisualAskUser (&plotFlag); + return true; +} + +bool pmReadoutCombineVisualCleanup(void) { + + if (!pmVisualIsVisual()) return true; + if (kapa == -1) return true; + + KapaClose(kapa); + return true; +} + +# else + +bool pmReadoutCombineVisualInit(void) { return true; } +bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean) { return true; } +bool pmReadoutCombineVisualCleanup(void) { return true; } + +# endif Index: branches/pap/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- branches/pap/psModules/src/imcombine/pmReadoutCombine.h (revision 23948) +++ branches/pap/psModules/src/imcombine/pmReadoutCombine.h (revision 25027) @@ -47,4 +47,8 @@ ); +bool pmReadoutCombineVisualInit(void); +bool pmReadoutCombineVisualPixels(psVector *pixels, psVector *mask, float mean); +bool pmReadoutCombineVisualCleanup(void); + /// @} #endif Index: branches/pap/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtraction.c (revision 25027) @@ -196,49 +196,12 @@ { psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image - int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy for (int y = -footprint; y <= footprint; y++) { - // Convolution with a delta function is just the value specified by the offset - memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes); + for (int x = -footprint; x <= footprint; x++) { + convolved->kernel[y][x] = image->kernel[y - v][x - u]; + } } return convolved; } -// Take a subset of an image -static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved - psRegion region, // Region of interest (with border) - bool threaded // Are we running threaded? - ) -{ - if (!image) { - return NULL; - } - // XXX if (threaded) { - // XXX psMutexLock(image); - // XXX } - psImage *subset = psImageSubset(image, region); // Subset image, to return - // XXX if (threaded) { - // XXX psMutexUnlock(image); - // XXX } - return subset; -} - -// Free the subset of an image -static inline void convolveSubsetFree(psImage *parent, // Parent image - psImage *child, // Child (subset) image - bool threaded // Are we running threaded? - ) -{ - if (!child || !parent) { - return; - } - // XXX if (threaded) { - // XXX psMutexLock(parent); - // XXX } - psFree(child); - // XXX if (threaded) { - // XXX psMutexUnlock(parent); - // XXX } - return; -} // Convolve an image using FFT @@ -256,13 +219,11 @@ region.y0 - size, region.y1 + size); // Add a border - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - - psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve - psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask + psImage *subImage = image ? psImageSubset(image, border) : NULL; // Subimage to convolve + psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution - convolveSubsetFree(image, subImage, threaded); - convolveSubsetFree(mask, subMask, threaded); + psFree(subImage); + psFree(subMask); // Now, we have to stick it in where it belongs @@ -300,9 +261,7 @@ region.y0 - size, region.y1 + size); // Add a border - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - - psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map - psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image - psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask + psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map + psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image + psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once @@ -310,7 +269,7 @@ psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys - convolveSubsetFree(variance, subVariance, threaded); - convolveSubsetFree(sys, subSys, threaded); - convolveSubsetFree(mask, subMask, threaded); + psFree(subVariance); + psFree(subSys); + psFree(subMask); // Now, we have to stick it in where it belongs @@ -420,15 +379,13 @@ if (box > 0) { int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - psRegion region = psRegionSet(colMin - box, colMax + box, rowMin - box, rowMax + box); // Region to convolve - psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve + psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad, -box, box, -box, box); // Convolved subtraction mask - convolveSubsetFree(subMask, image, threaded); + psFree(image); psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns"); @@ -653,7 +610,7 @@ float value = 0.0; // Value of convolved pixel int uMin = x - size, uMax = x + size; // Range for u - psF32 *xKernelData = xKernel->data.F32; // Kernel values + psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values psF32 *imageData = &image->kernel[y][uMin]; // Image values - for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) { + for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { value += *xKernelData * *imageData; } @@ -668,7 +625,7 @@ float value = 0.0; // Value of convolved pixel int vMin = y - size, vMax = y + size; // Range for v - psF32 *yKernelData = yKernel->data.F32; // Kernel values + psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! - for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) { + for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { value += *yKernelData * *imageData; } @@ -821,4 +778,10 @@ } psFree(mask); + + // XXX raise an error? + if (isnan(stats->sampleMean)) { + psFree(stats); + return -1; + } double mean, rms; // Mean and RMS of deviations Index: branches/pap/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (revision 25027) @@ -15,5 +15,7 @@ #include "pmSubtractionVisual.h" -//#define TESTING +// #define TESTING // TESTING output for debugging; may not work with threads! + +#define USE_VARIANCE // Include variance in equation? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -31,5 +33,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = image1->kernel[y][x] * image2->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sum += value; } } @@ -195,5 +201,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumC += value; } } @@ -218,5 +228,8 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; + double invNoise2 = 1.0; +#ifdef USE_VARIANCE + invNoise2 /= variance->kernel[y][x]; +#endif double value = input->kernel[y][x] * invNoise2; sumI += value; @@ -277,5 +290,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = target->kernel[y][x] * conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumTC += value; } } @@ -297,5 +314,8 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = target->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif sumIT += value * input->kernel[y][x]; sumT += value; @@ -366,5 +386,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumC += value; } } @@ -765,5 +789,23 @@ for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int iy = 0; iy < stamp->matrix1->numRows; iy++) { + for (int ix = 0; ix < stamp->matrix1->numCols; ix++) { + if (!isfinite(stamp->matrix1->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix1\n"); + } + } + } + for (int ix = 0; ix < stamp->vector1->n; ix++) { + if (!isfinite(stamp->vector1->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector1\n"); + } + } +#endif + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); @@ -774,7 +816,21 @@ } } + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector1\n"); + } + } +#endif + calculatePenalty(sumVector, kernels); #ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector1\n"); + } + } { psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); @@ -798,5 +854,5 @@ psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); psFree(sumMatrix); if (!luMatrix) { @@ -807,5 +863,14 @@ return NULL; } - kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation); + kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation); + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int ix = 0; ix < kernels->solution1->n; ix++) { + if (!isfinite(kernels->solution1->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector1\n"); + } + } +#endif psFree(sumVector); @@ -1062,4 +1127,23 @@ source = stamp->image1; convolutions = stamp->convolutions1; + + // Having convolved image1 and changed its normalisation, we need to renormalise the residual + // so that it is on the scale of image1. + psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm, + false); // Kernel image + if (!image) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel."); + return false; + } + double sumKernel = 0; // Sum of kernel, for normalising residual + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel + for (int y = 0; y < fullSize; y++) { + for (int x = 0; x < fullSize; x++) { + sumKernel += image->data.F32[y][x]; + } + } + psFree(image); + devNorm = 1.0 / sumKernel / numPixels; break; case PM_SUBTRACTION_MODE_2: Index: branches/pap/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtractionKernels.c (revision 25027) @@ -145,10 +145,8 @@ // Calculate moments - double sum = 0.0; // Sum of kernel component, for normalisation double moment = 0.0; // Moment, for penalty for (int v = -size, y = 0; v <= size; v++, y++) { for (int u = -size, x = 0; u <= size; u++, x++) { double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel - sum += value; moment += value * (PS_SQR(u) + PS_SQR(v)); } @@ -163,8 +161,8 @@ } } - sum = 1.0 / sum; + sum = 1.0 / sqrt(sum); psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); - moment *= sum; + moment *= PS_SQR(sum); } Index: branches/pap/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionMask.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtractionMask.c (revision 25027) @@ -6,4 +6,5 @@ #include +#include "pmErrorCodes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -78,5 +79,5 @@ } if (numBad > badFrac * numCols * numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, + psError(PM_ERR_SMALL_AREA, true, "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); Index: branches/pap/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionMatch.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtractionMatch.c (revision 25027) @@ -246,5 +246,5 @@ badFrac, mode); // Subtraction mask if (!subMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psError(psErrorCodeLast(), false, "Unable to generate subtraction mask."); psFree(kernels); psFree(regions); @@ -337,5 +337,5 @@ PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false); } - PS_ASSERT_INT_POSITIVE(iter, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); @@ -379,5 +379,5 @@ badFrac, subMode); if (!subMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psError(psErrorCodeLast(), false, "Unable to generate subtraction mask."); goto MATCH_ERROR; } @@ -432,10 +432,10 @@ } - if (sources) { + if (stampsName && strlen(stampsName) > 0) { + stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, + footprint, stampSpacing, subMode); + } else if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, footprint, stampSpacing, subMode); - } else if (stampsName && strlen(stampsName) > 0) { - stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, - footprint, stampSpacing, subMode); } @@ -830,12 +830,14 @@ psFree(mask); + // XXX raise an error here or not? + if (isnan(stats->robustMedian)) { + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median width ratio: %lf", stats->robustMedian); pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2); psFree(stats); - // XXX EAM : I think Paul left some test code in here. I've commented these lines out - // return PM_SUBTRACTION_MODE_2; - // exit(1); - return mode; } Index: branches/pap/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/pap/psModules/src/imcombine/pmSubtractionStamps.c (revision 23948) +++ branches/pap/psModules/src/imcombine/pmSubtractionStamps.c (revision 25027) @@ -78,29 +78,50 @@ } -// Is this position unmasked? -static bool checkStampMask(int x, int y, // Coordinates of stamp - int numCols, int numRows, // Size of image - const psImage *mask, // Mask - pmSubtractionMode mode, // Mode for subtraction - int footprint, // Size of stamp - int border // Size of border - ) -{ - // Check the footprint bounds - if (x < border || x >= numCols - border || y < border || y >= numRows - border) { - return false; - } - - if (!mask) { - return true; - } - - // Check the central pixel - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; - return false; - } - - return true; + +// Search a region for a suitable stamp +bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return + float *fluxStamp, // Flux of stamp, to return + const psImage *image1, const psImage *image2, // Images to search + float thresh1, float thresh2, // Thresholds for images + const psImage *subMask, // Subtraction mask + int xMin, int xMax, int yMin, int yMax, // Bounds of search + int numCols, int numRows, // Size of images + int border // Border around image + ) +{ + bool found = false; // Found a suitable stamp? + *fluxStamp = -INFINITY; // Flux of best stamp + + // Ensure we're not going to go outside the bounds of the image + xMin = PS_MAX(border, xMin); + xMax = PS_MIN(numCols - border - 1, xMax); + yMin = PS_MAX(border, yMin); + yMax = PS_MIN(numRows - border - 1, yMax); + + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + if ((image1 && image1->data.F32[y][x] < thresh1) || + (image2 && image2->data.F32[y][x] < thresh2)) { + continue; + } + + if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & + (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { + return false; + } + + // We take the MIN to attempt to avoid transients in both images + float flux = (image1 && image2) ? PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]) : + ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel + if (flux > *fluxStamp) { + *xStamp = x; + *yStamp = y; + *fluxStamp = flux; + found = true; + } + } + } + + return found; } @@ -314,7 +335,7 @@ numSearch++; - float xStamp = 0, yStamp = 0; // Coordinates of stamp - float fluxStamp = NAN; // Flux of stamp - bool goodStamp = false; // Found a good stamp? + int xStamp = 0, yStamp = 0; // Coordinates of stamp + float fluxStamp = -INFINITY; // Flux of stamp + bool goodStamp = false; // Found a good stamp? // A couple different ways of finding stamps: @@ -324,66 +345,35 @@ psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes - // Take stamp off the top of the (sorted) list - if (xList->n > 0) { - int index = xList->n - 1; // Index of new stamp - xStamp = xList->data.F32[index]; - yStamp = yList->data.F32[index]; - fluxStamp = fluxList->data.F32[index]; + // Take stamps off the top of the (sorted) list + for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { + int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre // Chop off the top of the list - xList->n = index; - yList->n = index; - fluxList->n = index; - - goodStamp = true; - int x = xStamp + 0.5, y = yStamp + 0.5; // Pixel location - if (image1 && image1->data.F32[y][x] < thresh1) { - goodStamp = false; - } else if (image2 && image2->data.F32[y][x] < thresh2) { - goodStamp = false; - } - } else { - psTrace("psModules.imcombine", 9, "No sources in subregion %d", i); + xList->n = j; + yList->n = j; + fluxList->n = j; + + // Fish around a bit to see if we can find a pixel that isn't masked + psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n", + i, xCentre, yCentre); + + // Search bounds + int search = footprint - size; // Search radius + int xMin = PS_MAX(border, xCentre - search); + int xMax = PS_MIN(numCols - border -1, xCentre + search); + int yMin = PS_MAX(border, yCentre - search); + int yMax = PS_MIN(numRows - border - 1, yCentre + search); + + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); } } else { // Use a simple method of automatically finding stamps --- take the highest pixel in the // subregion - fluxStamp = 0.0; psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest - const psImage *image = NULL; // Image for flux of stamp - switch (mode) { - case PM_SUBTRACTION_MODE_1: - image = image2 ? image2 : image1; - break; - case PM_SUBTRACTION_MODE_2: - image = image1 ? image1 : image2; - break; - case PM_SUBTRACTION_MODE_UNSURE: - case PM_SUBTRACTION_MODE_DUAL: - // If both are available, we'll take the MIN value below in the hope of rejecting - // transients in both - image = (image1 && image2) ? NULL : (image1 ? image1 : image2); - default: - psAbort("Unrecognised mode: %x", mode); - } - - for (int y = subRegion->y0; y <= subRegion->y1; y++) { - for (int x = subRegion->x0; x <= subRegion->x1; x++) { - if ((image1 && image1->data.F32[y][x] < thresh1) || - (image2 && image2->data.F32[y][x] < thresh2)) { - continue; - } - float value = image ? image->data.F32[y][x] : - PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]); // Value of image - if (value > fluxStamp && - checkStampMask(x, y, numCols, numRows, subMask, mode, footprint, border)) { - fluxStamp = image->data.F32[y][x]; - xStamp = x; - yStamp = y; - goodStamp = true; - } - } - } + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1, + numCols, numRows, border); } @@ -453,6 +443,4 @@ region, footprint, spacing); // Stamp list int numStamps = stamps->num; // Number of stamps - int numCols = image->numCols, numRows = image->numRows; // Size of image - int border = footprint + size; // Border size psString ds9name = NULL; // Filename for ds9 region file @@ -480,11 +468,4 @@ psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region", xPix, yPix); - pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta"); - continue; - } - if (!checkStampMask(xPix, yPix, numCols, numRows, subMask, mode, footprint, border)) { - // Not a good stamp - psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask", - xPix, yPix); pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red"); continue; Index: branches/pap/psModules/src/objects/Makefile.am =================================================================== --- branches/pap/psModules/src/objects/Makefile.am (revision 23948) +++ branches/pap/psModules/src/objects/Makefile.am (revision 25027) @@ -10,4 +10,5 @@ pmFootprintArraysMerge.c \ pmFootprintAssignPeaks.c \ + pmFootprintCullPeaks.c \ pmFootprintFind.c \ pmFootprintFindAtPoint.c \ @@ -38,4 +39,6 @@ pmSourceIO_PS1_CAL_0.c \ pmSourceIO_CMF_PS1_V1.c \ + pmSourceIO_CMF_PS1_V2.c \ + pmSourceIO_MatchedRefs.c \ pmSourcePlots.c \ pmSourcePlotPSFModel.c \ Index: branches/pap/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- branches/pap/psModules/src/objects/models/pmModel_PS1_V1.c (revision 23948) +++ branches/pap/psModules/src/objects/models/pmModel_PS1_V1.c (revision 25027) @@ -159,5 +159,5 @@ break; case PM_PAR_7: - params_min = 0.1; + params_min = -1.0; break; default: Index: branches/pap/psModules/src/objects/models/pmModel_SGAUSS.c =================================================================== --- branches/pap/psModules/src/objects/models/pmModel_SGAUSS.c (revision 23948) +++ branches/pap/psModules/src/objects/models/pmModel_SGAUSS.c (revision 25027) @@ -325,5 +325,8 @@ psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, contour, NULL, NULL, 0); + if (!psVectorStats (stats, contour, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } value = stats->sampleMedian; Index: branches/pap/psModules/src/objects/pmFootprint.h =================================================================== --- branches/pap/psModules/src/objects/pmFootprint.h (revision 23948) +++ branches/pap/psModules/src/objects/pmFootprint.h (revision 25027) @@ -51,12 +51,16 @@ void pmSetFootprintArrayIDsForImage(psImage *idImage, const psArray *footprints, // the footprints to insert - const bool relativeIDs); // show IDs starting at 0, not pmFootprint->id + const bool relativeIDs // show IDs starting at 0, not pmFootprint->id + ); psErrorCode pmFootprintsAssignPeaks(psArray *footprints, const psArray *peaks); -psErrorCode pmFootprintArrayCullPeaks(const psImage *img, const psImage *weight, psArray *footprints, - const float nsigma, const float threshold_min); -psErrorCode pmFootprintCullPeaks(const psImage *img, const psImage *weight, pmFootprint *fp, - const float nsigma, const float threshold_min); +psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + pmFootprint *fp, // Footprint containing mortal peaks + const float nsigma_delta, // how many sigma above local background a peak needs to be to survive + const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance + const float min_threshold // minimum permitted coll height + ); psArray *pmFootprintArrayToPeaks(const psArray *footprints); Index: branches/pap/psModules/src/objects/pmFootprintArraysMerge.c =================================================================== --- branches/pap/psModules/src/objects/pmFootprintArraysMerge.c (revision 23948) +++ branches/pap/psModules/src/objects/pmFootprintArraysMerge.c (revision 25027) @@ -28,18 +28,26 @@ */ psArray *pmFootprintArraysMerge(const psArray *footprints1, // one set of footprints - const psArray *footprints2, // the other set - const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2 - assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); - assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); + const psArray *footprints2, // the other set + const int includePeaks // which peaks to set? 0x1 => footprints1, 0x2 => 2 + ) +{ + if (!footprints1 && !footprints2) { + // No footprints in merged array + return psArrayAllocEmpty(0); + } - if (footprints1->n == 0 || footprints2->n == 0) { // nothing to do but put copies on merged - const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1; + assert(!footprints1 || footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); + assert(!footprints2 || footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); - psArray *merged = psArrayAllocEmpty(old->n); - for (int i = 0; i < old->n; i++) { - psArrayAdd(merged, 1, old->data[i]); - } - - return merged; + if (!footprints1 || footprints1->n == 0 || !footprints2 || footprints2->n == 0) { + // nothing to do but put copies on merged + const psArray *old = (!footprints1 || footprints1->n == 0) ? footprints2 : footprints1; + + psArray *merged = psArrayAllocEmpty(old->n); + for (int i = 0; i < old->n; i++) { + psArrayAdd(merged, 1, old->data[i]); + } + + return merged; } /* @@ -48,14 +56,14 @@ */ { - pmFootprint *fp1 = footprints1->data[0]; - pmFootprint *fp2 = footprints2->data[0]; - if (fp1->region.x0 != fp2->region.x0 || - fp1->region.x1 != fp2->region.x1 || - fp1->region.y0 != fp2->region.y0 || - fp1->region.y1 != fp2->region.y1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "The two pmFootprint arrays correspnond to different-sized regions"); - return NULL; - } + pmFootprint *fp1 = footprints1->data[0]; + pmFootprint *fp2 = footprints2->data[0]; + if (fp1->region.x0 != fp2->region.x0 || + fp1->region.x1 != fp2->region.x1 || + fp1->region.y0 != fp2->region.y0 || + fp1->region.y1 != fp2->region.y1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "The two pmFootprint arrays correspnond to different-sized regions"); + return NULL; + } } /* @@ -72,17 +80,17 @@ * Now assign the peaks appropriately. We could do this more efficiently * using idImage (which we just freed), but this is easy and probably fast enough - */ + */ if (includePeaks & 0x1) { - psArray *peaks = pmFootprintArrayToPeaks(footprints1); - pmFootprintsAssignPeaks(merged, peaks); - psFree(peaks); + psArray *peaks = pmFootprintArrayToPeaks(footprints1); + pmFootprintsAssignPeaks(merged, peaks); + psFree(peaks); } if (includePeaks & 0x2) { - psArray *peaks = pmFootprintArrayToPeaks(footprints2); - pmFootprintsAssignPeaks(merged, peaks); - psFree(peaks); + psArray *peaks = pmFootprintArrayToPeaks(footprints2); + pmFootprintsAssignPeaks(merged, peaks); + psFree(peaks); } - + return merged; } Index: branches/pap/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- branches/pap/psModules/src/objects/pmFootprintCullPeaks.c (revision 25027) +++ branches/pap/psModules/src/objects/pmFootprintCullPeaks.c (revision 25027) @@ -0,0 +1,265 @@ +/* @file pmFootprint.c + * low-level pmFootprint functions + * + * @author RHL, Princeton & IfA; EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-12-08 02:51:14 $ + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" + + /* + * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently + * isolated. More precisely, for each peak find the highest coll that you'd have to traverse + * to reach a still higher peak --- and if that coll's more than nsigma DN below your + * starting point, discard the peak. + */ + +# define IN_PEAK 1 +psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + pmFootprint *fp, // Footprint containing mortal peaks + const float nsigma_delta, // how many sigma above local background a peak needs to be to survive + const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance + const float min_threshold) { // minimum permitted coll height + assert (img != NULL); assert (img->type.type == PS_TYPE_F32); + assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); + assert (img->row0 == weight->row0 && img->col0 == weight->col0); + assert (fp != NULL); + + if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do + return PS_ERR_NONE; + } + + psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) + subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; + subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; + + psImage *subImg = psImageSubset((psImage *)img, subRegion); + psImage *subWt = psImageSubset((psImage *)weight, subRegion); + assert (subImg != NULL && subWt != NULL); + + psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); + + // We need a psArray of peaks brighter than the current peak. + // We reject peaks which either: + // 1) are below the local threshold + // 2) have a brighter peak within their threshold + + // allocate the full-sized array. if the final array is much smaller, we can realloc + // at that point. + psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n); + psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); + + // The brightest peak is always safe; go through other peaks trying to cull them + for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop + const pmPeak *peak = fp->peaks->data[i]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + // + // Find the level nsigma below the peak that must separate the peak + // from any of its friends + // + assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); + + // const float stdev = sqrt(subWt->data.F32[y][x]); + // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; + + const float stdev = sqrt(subWt->data.F32[y][x]); + const float flux = subImg->data.F32[y][x]; + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + // if flux is negative, careful with fStdev + + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold) || threshold < min_threshold) { + // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + continue; + } + + // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? + if (threshold > subImg->data.F32[y][x]) { + threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; + } + + // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint + // perhaps this should alloc a single ID image above and pass it in to be set. + + // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) + + // XXX this is not quite there yet: + // pmFootprint *peakFootprint = NULL; + // int area = subImg->numCols * subImg->numRows; + // if (area > 30000) { + + pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); + + // at this point brightPeaks only has the peaks brighter than the current + + // XXX need to supply the image here + // we set the IDs to either 1 (in peak) or 0 (not in peak) + pmSetFootprintID (idImg, peakFootprint, IN_PEAK); + psFree(peakFootprint); + + // Check if any of the previous (brighter) peaks are within the footprint of this peak + // If so, the current peak is bogus; drop it. + bool keep = true; + for (int j = 0; keep && (j < brightPeaks->n); j++) { + const pmPeak *peak2 = fp->peaks->data[j]; + int x2 = peak2->x - subImg->col0; + int y2 = peak2->y - subImg->row0; + if (idImg->data.S32[y2][x2] == IN_PEAK) + // There's a brighter peak within the footprint above threshold; so cull our initial peak + keep = false; + } + if (!keep) continue; + + psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + } + + psFree (fp->peaks); + fp->peaks = brightPeaks; + + psFree(idImg); + psFree(subImg); + psFree(subWt); + + return PS_ERR_NONE; +} + + + /* + * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently + * isolated. More precisely, for each peak find the highest coll that you'd have to traverse + * to reach a still higher peak --- and if that coll's more than nsigma DN below your + * starting point, discard the peak. + */ +psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + pmFootprint *fp, // Footprint containing mortal peaks + const float nsigma_delta, // how many sigma above local background a peak needs to be to survive + const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance + const float min_threshold) { // minimum permitted coll height + assert (img != NULL); assert (img->type.type == PS_TYPE_F32); + assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); + assert (img->row0 == weight->row0 && img->col0 == weight->col0); + assert (fp != NULL); + + if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do + return PS_ERR_NONE; + } + + psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) + subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; + subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; + const psImage *subImg = psImageSubset((psImage *)img, subRegion); + const psImage *subWt = psImageSubset((psImage *)weight, subRegion); + assert (subImg != NULL && subWt != NULL); + // + // We need a psArray of peaks brighter than the current peak. We'll fake this + // by reusing the fp->peaks but lying about n. + // + // We do this for efficiency (otherwise I'd need two peaks lists), and we are + // rather too chummy with psArray in consequence. But it works. + // + psArray *brightPeaks = psArrayAlloc(0); + psFree(brightPeaks->data); + brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks + // + // The brightest peak is always safe; go through other peaks trying to cull them + // + for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop + const pmPeak *peak = fp->peaks->data[i]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + // + // Find the level nsigma below the peak that must separate the peak + // from any of its friends + // + assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); + + // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations) + // add additional padding in quadrature: + + const float stdev = sqrt(subWt->data.F32[y][x]); + const float flux = subImg->data.F32[y][x]; + const float fStdev = stdev/flux; + const float stdev_pad = flux * hypot(fStdev, fPadding); + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold) || threshold < min_threshold) { +#if 1 // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + // XXX mark peak to be dropped + (void)psArrayRemoveIndex(fp->peaks, i); + i--; // we moved everything down one + continue; +#else +#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once + threshold = min_threshold; +#endif + } + + // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? + if (threshold > subImg->data.F32[y][x]) { + threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; + } + + // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint + // perhaps this should alloc a single ID image above and pass it in to be set. + + const int peak_id = 1; // the ID for the peak of interest + brightPeaks->n = i; // only stop at a peak brighter than we are + + // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) + + pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); + brightPeaks->n = 0; // don't double free + psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); + psFree(peakFootprint); + + // Check if any of the previous (brighter) peaks are within the footprint of this peak + // If so, the current peak is bogus; drop it. + int j; + for (j = 0; j < i; j++) { + const pmPeak *peak2 = fp->peaks->data[j]; + int x2 = peak2->x - subImg->col0; + int y2 = peak2->y - subImg->row0; + const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak + + if (peak2_id == peak_id) { // There's a brighter peak within the footprint above + ; // threshold; so cull our initial peak + (void)psArrayRemoveIndex(fp->peaks, i); + i--; // we moved everything down one + break; + } + } + if (j == i) { + j++; + } + + psFree(idImg); + } + + brightPeaks->n = 0; psFree(brightPeaks); + psFree((psImage *)subImg); + psFree((psImage *)subWt); + + return PS_ERR_NONE; +} + Index: branches/pap/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- branches/pap/psModules/src/objects/pmGrowthCurveGenerate.c (revision 23948) +++ branches/pap/psModules/src/objects/pmGrowthCurveGenerate.c (revision 25027) @@ -89,5 +89,8 @@ psVectorAppend (values, growth->fitMag); } - psVectorStats (stats, values, NULL, NULL, 0); + if (!psVectorStats (stats, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } psf->growth->fitMag = stats->sampleMedian; @@ -104,5 +107,8 @@ psVectorAppend (values, growth->apMag->data.F32[i]); } - psVectorStats (stats, values, NULL, NULL, 0); + if (!psVectorStats (stats, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } psf->growth->apMag->data.F32[i] = stats->sampleMedian; } Index: branches/pap/psModules/src/objects/pmMoments.c =================================================================== --- branches/pap/psModules/src/objects/pmMoments.c (revision 23948) +++ branches/pap/psModules/src/objects/pmMoments.c (revision 25027) @@ -51,4 +51,5 @@ tmp->Sky = 0.0; tmp->nPixels = 0; + tmp->SN = 0; psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); Index: branches/pap/psModules/src/objects/pmPSF.c =================================================================== --- branches/pap/psModules/src/objects/pmPSF.c (revision 23948) +++ branches/pap/psModules/src/objects/pmPSF.c (revision 25027) @@ -75,4 +75,6 @@ options->poissonErrorsParams = true; + options->chiFluxTrend = true; + return options; } Index: branches/pap/psModules/src/objects/pmPSF.h =================================================================== --- branches/pap/psModules/src/objects/pmPSF.h (revision 23948) +++ branches/pap/psModules/src/objects/pmPSF.h (revision 25027) @@ -82,4 +82,5 @@ bool poissonErrorsParams; ///< use poission errors for model parameter fitting float radius; + bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? } pmPSFOptions; Index: branches/pap/psModules/src/objects/pmPSFtry.c =================================================================== --- branches/pap/psModules/src/objects/pmPSFtry.c (revision 23948) +++ branches/pap/psModules/src/objects/pmPSFtry.c (revision 25027) @@ -43,12 +43,12 @@ for (int j = 0; j < trend->map->map->numRows; j++) { - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); - } - fprintf (stderr, "\t\t\t"); - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); - } - fprintf (stderr, "\n"); + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); + } + fprintf (stderr, "\t\t\t"); + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); + } + fprintf (stderr, "\n"); } return true; @@ -63,9 +63,9 @@ float Wt = 0.0; for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (!isfinite(map->map->data.F32[j][i])) continue; - Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i]; - Wt += map->error->data.F32[j][i]; - } + for (int i = 0; i < map->map->numCols; i++) { + if (!isfinite(map->map->data.F32[j][i])) continue; + Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i]; + Wt += map->error->data.F32[j][i]; + } } @@ -75,9 +75,9 @@ // XXX for now, we are just replacing bad pixels with the Mean for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (isfinite(map->map->data.F32[j][i]) && - (map->error->data.F32[j][i] > 0.0)) continue; - map->map->data.F32[j][i] = Mean; - } + for (int i = 0; i < map->map->numCols; i++) { + if (isfinite(map->map->data.F32[j][i]) && + (map->error->data.F32[j][i] > 0.0)) continue; + map->map->data.F32[j][i] = Mean; + } } return true; @@ -174,12 +174,12 @@ pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { + if (!source->moments) { psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - if (!source->moments->nPixels) { + continue; + } + if (!source->moments->nPixels) { psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } + continue; + } source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); @@ -211,5 +211,5 @@ if (Next == 0) { - psError(PS_ERR_UNKNOWN, true, "No sources with good extended fits from which to determine PSF."); + psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); psFree(psfTry); return NULL; @@ -282,5 +282,5 @@ if (Npsf == 0) { - psError(PS_ERR_UNKNOWN, true, "No sources with good PSF fits after model is built."); + psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); psFree(psfTry); return NULL; @@ -311,19 +311,22 @@ // linear clipped fit of chisq trend vs flux - bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, mask, 0xff, chisq, NULL, flux); - psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean - psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev - - psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", - psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); - - psFree(flux); - psFree(mask); - psFree(chisq); - - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); - psFree(psfTry); - return NULL; + if (options->chiFluxTrend) { + bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, + mask, 0xff, chisq, NULL, flux); + psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean + psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev + + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", + psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); + + psFree(flux); + psFree(mask); + psFree(chisq); + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); + psFree(psfTry); + return NULL; + } } @@ -600,7 +603,7 @@ for (int i = 0; i < sources->n; i++) { - // skip any masked sources (failed to fit one of the model steps or get a magnitude) - if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + pmSource *source = sources->data[i]; assert (source->modelEXT); // all unmasked sources should have modelEXT @@ -628,8 +631,8 @@ for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) { - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } + // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) + for (int i = 0; i < mask->n; i++) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { break; @@ -643,12 +646,12 @@ } if (entryMin == -1) { - psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params"); + psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params"); return false; } - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } + // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) + for (int i = 0; i < mask->n; i++) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); @@ -659,8 +662,8 @@ psf->trendNy = trend->map->map->numRows; - // copy mask back to srcMask - for (int i = 0; i < mask->n; i++) { - srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } + // copy mask back to srcMask + for (int i = 0; i < mask->n; i++) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } psFree (mask); @@ -686,33 +689,33 @@ for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; // XXX we are using the same stats structure on each pass: do we need to re-init it? bool status = true; - trend = psf->params->data[PM_PAR_E0]; + trend = psf->params->data[PM_PAR_E0]; status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); - pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); - - trend = psf->params->data[PM_PAR_E1]; + pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); + + trend = psf->params->data[PM_PAR_E1]; status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); - pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); - - trend = psf->params->data[PM_PAR_E2]; + pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); + + trend = psf->params->data[PM_PAR_E2]; status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); - pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); if (!status) { @@ -755,11 +758,11 @@ if (psf->fieldNx > psf->fieldNy) { Nx = scale; - float AR = psf->fieldNy / (float) psf->fieldNx; - Ny = (int) (Nx * AR + 0.5); + float AR = psf->fieldNy / (float) psf->fieldNx; + Ny = (int) (Nx * AR + 0.5); Ny = PS_MAX (1, Ny); } else { Ny = scale; - float AR = psf->fieldNx / (float) psf->fieldNy; - Nx = (int) (Ny * AR + 0.5); + float AR = psf->fieldNx / (float) psf->fieldNy; + Nx = (int) (Ny * AR + 0.5); Nx = PS_MAX (1, Nx); } @@ -809,41 +812,41 @@ if (scale == 1) { - x_fit = psMemIncrRefCounter (x); - y_fit = psMemIncrRefCounter (y); - x_tst = psMemIncrRefCounter (x); - y_tst = psMemIncrRefCounter (y); - e0obs_fit = psMemIncrRefCounter (e0obs); - e1obs_fit = psMemIncrRefCounter (e1obs); - e2obs_fit = psMemIncrRefCounter (e2obs); - e0obs_tst = psMemIncrRefCounter (e0obs); - e1obs_tst = psMemIncrRefCounter (e1obs); - e2obs_tst = psMemIncrRefCounter (e2obs); + x_fit = psMemIncrRefCounter (x); + y_fit = psMemIncrRefCounter (y); + x_tst = psMemIncrRefCounter (x); + y_tst = psMemIncrRefCounter (y); + e0obs_fit = psMemIncrRefCounter (e0obs); + e1obs_fit = psMemIncrRefCounter (e1obs); + e2obs_fit = psMemIncrRefCounter (e2obs); + e0obs_tst = psMemIncrRefCounter (e0obs); + e1obs_tst = psMemIncrRefCounter (e1obs); + e2obs_tst = psMemIncrRefCounter (e2obs); } else { - x_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - x_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - for (int i = 0; i < e0obs_fit->n; i++) { - // e0obs->n == 8 or 9: - // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6 - // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7 - x_fit->data.F32[i] = x->data.F32[2*i+0]; - x_tst->data.F32[i] = x->data.F32[2*i+1]; - y_fit->data.F32[i] = y->data.F32[2*i+0]; - y_tst->data.F32[i] = y->data.F32[2*i+1]; - - e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; - e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; - e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0]; - e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1]; - e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0]; - e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1]; - } + x_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + y_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + x_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + y_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + for (int i = 0; i < e0obs_fit->n; i++) { + // e0obs->n == 8 or 9: + // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6 + // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7 + x_fit->data.F32[i] = x->data.F32[2*i+0]; + x_tst->data.F32[i] = x->data.F32[2*i+1]; + y_fit->data.F32[i] = y->data.F32[2*i+0]; + y_tst->data.F32[i] = y->data.F32[2*i+1]; + + e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; + e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; + e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0]; + e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1]; + e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0]; + e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1]; + } } @@ -852,5 +855,5 @@ // copy mask values to fitMask as a starting point for (int i = 0; i < fitMask->n; i++) { - fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; } @@ -859,42 +862,42 @@ for (int i = 0; i < nIter; i++) { // XXX we are using the same stats structure on each pass: do we need to re-init it? - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - - trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); - - trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); - - trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + + trend = psf->params->data[PM_PAR_E0]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); + + trend = psf->params->data[PM_PAR_E1]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); + + trend = psf->params->data[PM_PAR_E2]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); } psf->psfTrendStats->clipIter = nIter; // restore default setting @@ -941,5 +944,5 @@ // XXX copy fitMask values back to mask for (int i = 0; i < fitMask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; } psFree (fitMask); @@ -958,13 +961,22 @@ float dEsquare = 0.0; psStatsInit (stats); - psVectorStats (stats, e0res, NULL, mask, maskValue); + if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); psStatsInit (stats); - psVectorStats (stats, e1res, NULL, mask, maskValue); + if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); psStatsInit (stats); - psVectorStats (stats, e2res, NULL, mask, maskValue); + if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); @@ -998,5 +1010,5 @@ for (int i = 0; i < nBin; i++) { int j; - int nValid = 0; + int nValid = 0; for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { int N = index->data.U32[n]; @@ -1006,7 +1018,7 @@ mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N]; - if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; - } - if (nValid < 3) continue; + if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; + } + if (nValid < 3) continue; dE0subset->n = j; @@ -1018,13 +1030,20 @@ float dEsquare = 0.0; psStatsInit (statsS); - psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff); + if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) { + } dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); psStatsInit (statsS); - psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff); + if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); psStatsInit (statsS); - psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff); + if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); Index: branches/pap/psModules/src/objects/pmPeaks.c =================================================================== --- branches/pap/psModules/src/objects/pmPeaks.c (revision 23948) +++ branches/pap/psModules/src/objects/pmPeaks.c (revision 25027) @@ -147,5 +147,5 @@ tmp->y = y; tmp->value = value; - tmp->flux = 0; + tmp->flux = value; tmp->SN = 0; tmp->xf = x; Index: branches/pap/psModules/src/objects/pmSource.c =================================================================== --- branches/pap/psModules/src/objects/pmSource.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSource.c (revision 25027) @@ -242,5 +242,4 @@ extend |= (mySource->maskView == NULL); - // extend = true; if (extend) { // re-create the subimage @@ -250,5 +249,5 @@ mySource->pixels = psImageSubset(readout->image, newRegion); - mySource->variance = psImageSubset(readout->variance, newRegion); + mySource->variance = psImageSubset(readout->variance, newRegion); mySource->maskView = psImageSubset(readout->mask, newRegion); mySource->region = newRegion; @@ -291,7 +290,7 @@ bool status = true; // Status of MD lookup - float PSF_CLUMP_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_SN_LIM"); + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); if (!status) { - PSF_CLUMP_SN_LIM = 0; + PSF_SN_LIM = 0; } float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); @@ -342,5 +341,5 @@ } - if (src->moments->SN < PSF_CLUMP_SN_LIM) { + if (src->moments->SN < PSF_SN_LIM) { psTrace("psModules.objects", 10, "Rejecting source from clump because of low S/N (%f)\n", src->moments->SN); @@ -450,5 +449,5 @@ if (tmpSrc->moments == NULL) continue; - if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM) + if (tmpSrc->moments->SN < PSF_SN_LIM) continue; @@ -479,9 +478,15 @@ stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - psVectorStats (stats, tmpSx, NULL, NULL, 0); + if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats"); + return (emptyClump); + } psfClump.X = stats->clippedMean; psfClump.dX = stats->clippedStdev; - psVectorStats (stats, tmpSy, NULL, NULL, 0); + if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats"); + return (emptyClump); + } psfClump.Y = stats->clippedMean; psfClump.dY = stats->clippedStdev; @@ -528,10 +533,9 @@ bool status; float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM"); - if (!status) - PSF_SN_LIM = 20.0; + if (!status) PSF_SN_LIM = 20.0; float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - if (!status) - PSF_CLUMP_NSIGMA = 1.5; - float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + if (!status) PSF_CLUMP_NSIGMA = 1.5; + + // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; @@ -543,12 +547,12 @@ pmSource *source = (pmSource *) sources->data[i]; - // psf clumps are found for image subregions: - // skip sources not in this region + // psf clumps are found for image subregions: + // skip sources not in this region if (source->peak->x < region->x0) continue; if (source->peak->x >= region->x1) continue; if (source->peak->y < region->y0) continue; - if (source->peak->y >= region->y1) continue; - - // should be set by pmSourceAlloc + if (source->peak->y >= region->y1) continue; + + // should be set by pmSourceAlloc psAssert (source->type == PM_SOURCE_TYPE_UNKNOWN, "source type was not init-ed?"); @@ -556,5 +560,5 @@ if (!source->moments) { source->type = PM_SOURCE_TYPE_STAR; - psAssert (source->mode & noMoments, "why is this source missing moments?"); + psAssert (source->mode & noMoments, "why is this source missing moments?"); Nstar++; continue; @@ -576,6 +580,7 @@ source->type = PM_SOURCE_TYPE_STAR; source->mode |= PM_SOURCE_MODE_SATSTAR; - // recalculate moments here with larger box? - pmSourceMoments (source, INNER_RADIUS); + // why do we recalculate moments here? + // we already attempt to do this in psphotSourceStats + // pmSourceMoments (source, INNER_RADIUS); Nsatstar ++; continue; @@ -590,36 +595,38 @@ } - // likely defect (too small to be stellar) (push out to 3 sigma) - // low S/N objects which are small are probably stellar - // only set candidate defects if - // XXX these limits are quite arbitrary - if ((sigX < 0.05) || (sigY < 0.05)) { - source->type = PM_SOURCE_TYPE_DEFECT; - source->mode |= PM_SOURCE_MODE_DEFECT; - Ncr ++; - continue; - } - - // likely unsaturated extended source (too large to be stellar) - if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - source->type = PM_SOURCE_TYPE_EXTENDED; - Next ++; - continue; - } - - // the rest are probable stellar objects - starsn_moments->data.F32[starsn_moments->n] = source->moments->SN; - starsn_moments->n ++; - starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN; - starsn_peaks->n ++; - Nstar ++; - - // PSF star (within 1.5 sigma of clump center, S/N > limit) - psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); - if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) { - source->type = PM_SOURCE_TYPE_STAR; - source->mode |= PM_SOURCE_MODE_PSFSTAR; - Npsf ++; - continue; + // The following determinations require the use of moments + if (!(source->mode & noMoments)) { + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // XXX these limits are quite arbitrary + if (sigX < 0.05 || sigY < 0.05) { + source->type = PM_SOURCE_TYPE_DEFECT; + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // likely unsaturated extended source (too large to be stellar) + if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) { + source->type = PM_SOURCE_TYPE_EXTENDED; + Next ++; + continue; + } + + // the rest are probable stellar objects + starsn_moments->data.F32[starsn_moments->n] = source->moments->SN; + starsn_moments->n ++; + starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN; + starsn_peaks->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center, S/N > limit) + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) { + source->type = PM_SOURCE_TYPE_STAR; + source->mode |= PM_SOURCE_MODE_PSFSTAR; + Npsf ++; + continue; + } } @@ -636,6 +643,8 @@ if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) { - // Don't care about this error - psErrorClear(); + psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); + psFree (stats); + psFree (starsn_peaks); + return false; } psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats->min, stats->max); @@ -648,6 +657,8 @@ stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) { - // Don't care about this error - psErrorClear(); + psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); + psFree (stats); + psFree (starsn_peaks); + return false; } psLogMsg ("psModules.objects", 3, "SN range (peaks) : %f - %f (%ld)\n", Index: branches/pap/psModules/src/objects/pmSource.h =================================================================== --- branches/pap/psModules/src/objects/pmSource.h (revision 23948) +++ branches/pap/psModules/src/objects/pmSource.h (revision 25027) @@ -213,6 +213,8 @@ */ bool pmSourceMoments( - pmSource *source, ///< The input pmSource for which moments will be computed - float radius ///< Use a circle of pixels around the peak + pmSource *source, ///< The input pmSource for which moments will be computed + float radius, ///< Use a circle of pixels around the peak + float sigma, ///< size of Gaussian window function (<= 0.0 -> skip window) + float minSN ///< minimum pixel significance ); Index: branches/pap/psModules/src/objects/pmSourceContour.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceContour.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceContour.c (revision 25027) @@ -273,5 +273,5 @@ x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); if (x0 == x1) { - // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); goto pt1; } @@ -282,5 +282,5 @@ x1 = findContourNeg (image, threshold, xR, yR, RIGHT); } - // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); xVec->data.F32[Npt + 0] = image->col0 + x0; @@ -288,4 +288,5 @@ yVec->data.F32[Npt + 0] = image->row0 + yR; yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; @@ -307,5 +308,5 @@ x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); if (x0 == x1) { - // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); goto pt2; } @@ -322,4 +323,5 @@ yVec->data.F32[Npt + 0] = image->row0 + yR; yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; @@ -334,5 +336,4 @@ yVec->n = Npt; - // fprintf (stderr, "done\n"); psArray *tmpArray = psArrayAlloc(2); Index: branches/pap/psModules/src/objects/pmSourceFitModel.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceFitModel.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceFitModel.c (revision 25027) @@ -172,10 +172,10 @@ fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc); for (int i = 0; i < dparams->n; i++) { - if (psTraceGetLevel("psModules.objects") >= 4) { - fprintf (stderr, "%f ", params->data.F32[i]); - } if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); + if (psTraceGetLevel("psModules.objects") >= 4) { + fprintf (stderr, "%f +/- %f\n", params->data.F32[i], dparams->data.F32[i]); + } } psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); Index: branches/pap/psModules/src/objects/pmSourceGroup.h =================================================================== --- branches/pap/psModules/src/objects/pmSourceGroup.h (revision 25027) +++ branches/pap/psModules/src/objects/pmSourceGroup.h (revision 25027) @@ -0,0 +1,44 @@ +/* @file pmSourceGroup.h + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# ifndef PM_SOURCE_GROUP_H +# define PM_SOURCE_GROUP_H + +#include +#include "pmPeaks.h" +#include "pmModel.h" +#include "pmMoments.h" +#include "pmSourceExtendedPars.h" + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +#include + +/** pmSourceGroup data structure + * + * A source group is a connected set of source measurements with a + * common connection. Each source in the group represents the + * detection of some astronomical object on a single each. The group + * represents the related collection of measurements. The fits are + * coupled in some way. For example, they may all have the same + * position, but independent fluxes. Or, they may have a common set + * of positions and shape parameters. Or the position in each image + * may be related by a function. + * + * XXX is thre any info that is neaded for each source group beyond that carried by the sources (besides ID)? + */ +struct pmSource { + int seq; ///< ID for output (generated on write OR set on read) + psArray *sources; +} pmSourceGroup; + + +/// @} +# endif /* PM_SOURCE_GROUP_H */ Index: branches/pap/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceIO.c (revision 25027) @@ -496,4 +496,7 @@ status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); + } if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { @@ -506,4 +509,7 @@ status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); + } } if (xfitname) { @@ -517,4 +523,7 @@ status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); + } } if (!status) { @@ -562,6 +571,8 @@ // not needed if only one chip - if (file->fpa->chips->n == 1) return true; - + if (file->fpa->chips->n == 1) { + pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); + return true; + } // find the FPA phu @@ -667,4 +678,5 @@ psFree (outhead); + pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); return true; } @@ -679,4 +691,6 @@ pmFPA *fpa = file->fpa; + + pmSourceIO_ReadMatchedRefs (file->fits, fpa, config); if (view->chip == -1) { @@ -941,4 +955,7 @@ sources = pmSourcesRead_CMF_PS1_V1 (file->fits, hdu->header); } + if (!strcmp (exttype, "PS1_V2")) { + sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); + } } @@ -1053,2 +1070,3 @@ } + Index: branches/pap/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO.h (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceIO.h (revision 25027) @@ -39,4 +39,8 @@ bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname); + bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); @@ -48,4 +52,5 @@ psArray *pmSourcesRead_PS1_CAL_0 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); bool pmSourcesWritePSFs (psArray *sources, char *filename); @@ -75,4 +80,7 @@ bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos); +bool pmSourceIO_WriteMatchedRefs (psFits *fits, pmFPA *fpa, pmConfig *config); +bool pmSourceIO_ReadMatchedRefs (psFits *fits, pmFPA *fpa, const pmConfig *config); + /// @} # endif /* PM_SOURCE_IO_H */ Index: branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25027) +++ branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25027) @@ -0,0 +1,625 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.h" + +// panstarrs-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, + psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + psF32 errMag, chisq, apRadius; + psS32 nPix, nDOF; + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + bool status1 = false; + bool status2 = false; + float magOffset = NAN; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + + // if the sequence is defined, write these in seq order; otherwise + // write them in S/N order: + if (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state + // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // If source->seq is -1, source was generated in this analysis. If source->seq is + // not -1, source was read from elsewhere: in the latter case, preserve the source + // ID. source.seq is used instead of source.id since the latter is a const + // generated on Alloc, and would thus be wrong for read in sources. + if (source->seq == -1) { + source->seq = i; + } + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + } else { + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + } + chisq = model->chisq; + nDOF = model->nDOF; + nPix = model->nPix; + apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + chisq = NAN; + nDOF = 0; + nPix = 0; + apRadius = NAN; + errMag = NAN; + } + + float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); + + // distinguish moments measure from window vs S/N > XX ?? + float mxx = source->moments ? source->moments->Mxx : NAN; + float mxy = source->moments ? source->moments->Mxy : NAN; + float myy = source->moments ? source->moments->Myy : NAN; + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + // XXX not sure how to get this : need to load Nimages with weight? + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + + psArrayAdd (table, 100, row); + psFree (row); + + // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not + // subtracted + + // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image + // edge; 3) any pixels in the 3x3 peak region are masked; + } + + psMetadata *header = psMetadataCopy(NULL, tableHeader); + pmSourceMasksHeader(header); + + if (table->n == 0) { + psFitsWriteBlank(fits, header, extname); + psFree(table); + psFree(header); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + assert (modelType > -1); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entriesa + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN; + + pmPSF_AxesToModel (PAR, axes); + + float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + source->peak->flux = peakFlux; + source->peak->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + source->moments = pmMomentsAlloc (); + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + + bool status; + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + assert (radialBinsLower->n == radialBinsUpper->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + // skip source if it is not a ext sourc + // XXX we have two places that extended source parameters are measured: + // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, + // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) + // should we require both? + + pmSource *source = sources->data[i]; + + // skip sources without measurements + if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; + + // XXX I need to split the extended models from the extended aperture measurements + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + if (doPetrosian) { + pmSourcePetrosianValues *petrosian = source->extpars->petrosian; + if (petrosian) { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + } + } + + // Kron measurements + if (doKron) { + pmSourceKronValues *kron = source->extpars->kron; + if (kron) { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); + } + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { + pmSourceIsophotalValues *isophot = source->extpars->isophot; + if (isophot) { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); + } + } + + // Flux Annuli + if (doAnnuli) { + pmSourceAnnuli *annuli = source->extpars->annuli; + if (annuli) { + psVector *fluxVal = annuli->flux; + psVector *fluxErr = annuli->fluxErr; + psVector *fluxVar = annuli->fluxVar; + + for (int j = 0; j < fluxVal->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); + } + } else { + for (int j = 0; j < radialBinsLower->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); + } + } + } + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + + return true; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} Index: branches/pap/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 25027) +++ branches/pap/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 25027) @@ -0,0 +1,184 @@ +/** @file pmSourceIO_MatchedRefs.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.h" +#include "pmAstrometryObjects.h" +#include "pmAstrometryWCS.h" + +bool pmSourceIO_WriteMatchedRefs (psFits *fits, pmFPA *fpa, pmConfig *config) { + + bool status = true; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + + // first, check if there are any matches to be written + + // determine the output table format + // XXX move this elsewhere? (psastro recipe? filerules?) + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSASTRO"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSASTRO in config data"); + return false; + } + + bool REFS_OUTPUT = psMetadataLookupBool(&status, recipe, "PSASTRO.SAVE.REFMATCH"); + if (!REFS_OUTPUT) return true; + + psArray *table = psMetadataLookupPtr (&status, fpa->analysis, "MATCHED_REFS"); + psMemIncrRefCounter (table); + + if (!table) { + table = psArrayAllocEmpty (0x1000); + pmFPAview *view = pmFPAviewAlloc (0); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) continue; + + char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) continue; + + // process each of the readouts + // XXX there can only be one readout per chip, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) continue; + + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) continue; + + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) continue; + psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); + +# if (0) + // XXX test + FILE *outfile = fopen ("refstars.dat", "w"); + assert (outfile); + for (int nn = 0; nn < refstars->n; nn++) { + pmAstromObj *ref = refstars->data[nn]; + fprintf (outfile, "%lf %lf\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD); + } + fclose (outfile); +# endif + + psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (matches == NULL) continue; + + for (int i = 0; i < matches->n; i++) { + pmAstromMatch *match = matches->data[i]; + + pmAstromObj *raw = rawstars->data[match->raw]; + pmAstromObj *ref = refstars->data[match->ref]; + + psMetadata *row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "RA", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); + psMetadataAdd (row, PS_LIST_TAIL, "DEC", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color", ref->Color); + psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); + // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined. + + psArrayAdd (table, 100, row); + psFree (row); + } + } + } + } + psFree (view); + + if (table->n == 0) { + psFree(table); + return true; + } + } + + if (!psFitsWriteTable(fits, NULL, table, "MATCHED_REFS")) { + psError(PS_ERR_IO, false, "writing MATCHED_REFS\n"); + psFree(table); + return false; + } + + psFree(table); + return true; +} + +bool pmSourceIO_ReadMatchedRefs (psFits *fits, pmFPA *fpa, const pmConfig *config) { + + bool status = true; + + // check if we've already read (attempted to read) REFMATCH + bool readMatchedRefs = psMetadataLookupBool (&status, fpa->analysis, "READ.REFMATCH"); + if (readMatchedRefs) return true; + + // try find the MATCHED_REFS extension. if non-existent, note that we tried, and move on. + if (!psFitsMoveExtName (fits, "MATCHED_REFS")) { + psMetadataAddBool (fpa->analysis, PS_LIST_TAIL, "READ.REFMATCH", PS_META_REPLACE, "attempted to read MATCHED_REFS", true); + return true; + } + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numRows = psFitsTableSize(fits); // Number of sources in table + psArray *rows = psArrayAlloc(numRows); // Array of sources, to return + + // first, check if there are any matches to be written + for (int i = 0; i < numRows; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + rows->data[i] = row; + } + + psMetadataAddArray (fpa->analysis, PS_LIST_TAIL, "MATCHED_REFS", PS_META_REPLACE, "MATCHED_REFS", rows); + psFree (rows); + + return true; +} Index: branches/pap/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_RAW.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceIO_RAW.c (revision 25027) @@ -254,9 +254,12 @@ } + fprintf (f, "# %5s %5s %8s %7s %7s %6s %6s %10s %7s %7s %7s %4s %4s %5s\n", + "x", "y", "peak", "Mx", "My", "Mxx", "Myy", "Sum", "Peak", "Sky", "SN", "nPix", "type", "mode"); + for (i = 0; i < sources->n; i++) { source = sources->data[i]; if (source->moments == NULL) continue; - fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", + fprintf (f, "%5d %5d %8.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", source->peak->x, source->peak->y, source->peak->value, source->moments->Mx, source->moments->My, Index: branches/pap/psModules/src/objects/pmSourceMasks.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceMasks.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceMasks.c (revision 25027) @@ -44,4 +44,7 @@ ADD_MASK(header, RADIAL_FLUX , "Calculated radial flux measurements"); ADD_MASK(header, SIZE_SKIPPED , "Could not be determine size"); + ADD_MASK(header, ON_SPIKE , "Source lands in bright star spike"); + ADD_MASK(header, ON_GHOST , "Source lands in bright star ghost / glint"); + ADD_MASK(header, OFF_CHIP , "Source centroid lands off chip edge"); return true; Index: branches/pap/psModules/src/objects/pmSourceMasks.h =================================================================== --- branches/pap/psModules/src/objects/pmSourceMasks.h (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceMasks.h (revision 25027) @@ -34,4 +34,7 @@ PM_SOURCE_MODE_RADIAL_FLUX = 0x08000000, ///< radial flux measurements calculated PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip } pmSourceMode; Index: branches/pap/psModules/src/objects/pmSourceMatch.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceMatch.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceMatch.c (revision 25027) @@ -45,4 +45,5 @@ psVector **x, psVector **y, // Coordinate vectors to return psVector **mag, psVector **magErr, // Magnitude and error vectors to return + psVector **indices, // Indices for sources const psArray *sources // Input sources ) @@ -58,13 +59,13 @@ *mag = psVectorRecycle(*mag, numSources, PS_TYPE_F32); *magErr = psVectorRecycle(*magErr, numSources, PS_TYPE_F32); + *indices = psVectorRecycle(*indices, numSources, PS_TYPE_S32); float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources long num = 0; // Number of valid sources for (long i = 0; i < numSources; i++) { pmSource *source = sources->data[i]; // Source of interest - if (!source) continue; - if (source->mode & SOURCE_MASK) continue; - if (!isfinite(source->psfMag)) continue; - if (!isfinite(source->errMag)) continue; - if (source->psfMag > SOURCE_FAINTEST) continue; + if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || + !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) { + continue; + } float xSrc, ySrc; // Coordinates of source @@ -78,5 +79,6 @@ (*y)->data.F32[num] = ySrc; (*mag)->data.F32[num] = source->psfMag; - (*magErr)->data.F32[num] = source->errMag ; + (*magErr)->data.F32[num] = source->errMag; + (*indices)->data.S32[num] = i; num++; } @@ -85,4 +87,5 @@ (*mag)->n = num; (*magErr)->n = num; + (*indices)->n = num; if (*bounds) { @@ -175,6 +178,7 @@ psVector *xImage = NULL, *yImage = NULL; // Coordinates of sources psVector *magImage = NULL, *magErrImage = NULL; // Magnitude and mag - - int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, + psVector *indices = NULL; // Indices for sources + + int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, &indices, sources); // Number of sources @@ -187,5 +191,6 @@ for (int j = 0; j < numSources; j++) { pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); matches->data[j] = match; } @@ -206,5 +211,6 @@ for (int j = 0, k = numMaster; j < numSources; j++, k++) { pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); matches->data[k] = match; } @@ -231,10 +237,12 @@ // Record the match pmSourceMatch *match = matches->data[index]; // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); numMatch++; } else { // Add to the master list pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); xMaster->data.F32[numMaster] = xImage->data.F32[j]; yMaster->data.F32[numMaster] = yImage->data.F32[j]; @@ -469,4 +477,10 @@ if (!psVectorStats(stats, trans, NULL, badImage, 0xFF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies."); + psFree(stats); + return -1; + } + // XXX handle this case better: + if (isnan(stats->clippedMean)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies."); psFree(stats); Index: branches/pap/psModules/src/objects/pmSourceMoments.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceMoments.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceMoments.c (revision 25027) @@ -57,5 +57,5 @@ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) -bool pmSourceMoments(pmSource *source, psF32 radius) +bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -84,4 +84,6 @@ psF32 Y1 = 0.0; psF32 R2 = PS_SQR(radius); + psF32 minSN2 = PS_SQR(minSN); + psF32 rsigma2 = 0.5 / PS_SQR(sigma); // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -97,4 +99,7 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yPeak; + if (fabs(yDiff) > radius) continue; psF32 *vPix = source->pixels->data.F32[row]; @@ -113,5 +118,5 @@ psF32 xDiff = col - xPeak; - psF32 yDiff = row - yPeak; + if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) @@ -121,15 +126,33 @@ psF32 wDiff = *vWgt; - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; + // skip pixels below specified significance level + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + if (pDiff < 0) continue; + + if (sigma > 0.0) { + // apply a pseudo-gaussian weight + + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); + psF32 weight = 1.0 / t; + + // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); + + wDiff *= weight; + pDiff *= weight; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + // fprintf (stderr, " : %6.1f %6.1f %8.1f %8.1f\n", X1, Y1, Sum, Var); peakPixel = PS_MAX (*vPix, peakPixel); @@ -188,4 +211,7 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radius) continue; + psF32 *vPix = source->pixels->data.F32[row]; psF32 *vWgt = source->variance->data.F32[row]; @@ -203,5 +229,5 @@ psF32 xDiff = col - xCM; - psF32 yDiff = row - yCM; + if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) @@ -215,6 +241,21 @@ // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; + if (PS_SQR(pDiff) < minSN2*wDiff) continue; if (pDiff < 0) continue; + + if (sigma > 0.0) { + // apply a pseudo-gaussian weight + + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); + psF32 weight = 1.0 / t; + + // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); + + wDiff *= weight; + pDiff *= weight; + } Sum += pDiff; Index: branches/pap/psModules/src/objects/pmSourceUtils.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceUtils.c (revision 23948) +++ branches/pap/psModules/src/objects/pmSourceUtils.c (revision 25027) @@ -94,24 +94,35 @@ source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE); - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + float xReadout, yReadout; - // XXX fix the binning : currently not selected from concepts - // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y - // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y - int xBin = 1; - int yBin = 1; + // if we have information about the chip & cell, adjust the coordinates chip->cell->readout + // otherwise, assume 0,0 offset and 1,1 parity + if (cell) { + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - // Position on the cell - float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); - float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); + // XXX fix the binning : currently not selected from concepts + // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + int xBin = 1; + int yBin = 1; - // Position on the readout - // float xReadout = CELL_TO_READOUT(xCell, x0Readout); - // float yReadout = CELL_TO_READOUT(yCell, y0Readout); + // Position on the cell + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); + + // Position on the readout + // float xReadout = CELL_TO_READOUT(xCell, x0Readout); + // float yReadout = CELL_TO_READOUT(yCell, y0Readout); + xReadout = xCell; + yReadout = yCell; + } else { + xReadout = xChip; + yReadout = yChip; + } - pmSourceDefinePixels (source, readout, xCell, yCell, radius); + pmSourceDefinePixels (source, readout, xReadout, yReadout, radius); return (source); Index: branches/pap/psModules/src/psmodules.h =================================================================== --- branches/pap/psModules/src/psmodules.h (revision 23948) +++ branches/pap/psModules/src/psmodules.h (revision 25027) @@ -47,4 +47,5 @@ #include #include +#include #include #include @@ -76,4 +77,5 @@ #include #include +#include // the following headers are from psModule:astrom Index: branches/pap/psModules/test/camera/tap_pmFPAMaskW.c =================================================================== --- branches/pap/psModules/test/camera/tap_pmFPAMaskW.c (revision 23948) +++ branches/pap/psModules/test/camera/tap_pmFPAMaskW.c (revision 25027) @@ -356,5 +356,5 @@ // ---------------------------------------------------------------------- // pmReadoutSetVariance() tests: NULL inputs - // bool pmReadoutSetVariance(pmReadout *readout, bool poisson) + // bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson) if (1) { psMemId id = psMemGetId(); @@ -367,6 +367,6 @@ // Set readout == NULL, ensure pmReadoutSetVariance() returnes FALSE, with no seg faults, memory leaks - rc = pmReadoutSetVariance(NULL, false); - ok(!rc, "pmReadoutSetVariance(NULL, false) returned FALSE with null pmReadout input"); + rc = pmReadoutSetVariance(NULL, NULL, false); + ok(!rc, "pmReadoutSetVariance(NULL, NULL, false) returned FALSE with null pmReadout input"); @@ -375,5 +375,7 @@ rc|= psMetadataAddF32(readout->parent->concepts, PS_LIST_HEAD, "CELL.READNOISE", PS_META_REPLACE, NULL, CELL_READNOISE); ok(rc, "Set GAIN and READNOISE in cell->concepts successfully"); - +/* + * Getting the section below to run requires generating a noiseMap + * // Call pmReadoutSetVariance() and then verify that the mask data was set correctly rc = pmReadoutSetVariance(readout, false); @@ -419,4 +421,5 @@ ok(!errorFlag, "pmReadoutSetWeight() set the weight values correctly (Poisson)"); +*/ psFree(fpa); psFree(camera); Index: branches/pap/psModules/test/concepts/tap_pmConcepts.c =================================================================== --- branches/pap/psModules/test/concepts/tap_pmConcepts.c (revision 23948) +++ branches/pap/psModules/test/concepts/tap_pmConcepts.c (revision 25027) @@ -77,4 +77,22 @@ } +psMetadataItem *dummyConceptCopier( + const psMetadataItem *source, + const psMetadataItem *target, + const psMetadata *cameraFormat, + const pmFPA *fpa, + const pmChip *chip, + const pmCell *cell) +{ + if (target == NULL || + source == PM_CONCEPT_SOURCE_NONE || + cameraFormat == NULL || + fpa == NULL || + chip == NULL || + cell == NULL) { + printf("dummyConceptCopier() args are NULL\n"); + } + return(NULL); +} int main(int argc, char* argv[]) @@ -92,5 +110,5 @@ psMetadataItem *blank = psMetadataItemAlloc("myItem1", PS_DATA_BOOL, "I am a boolean", true); pmConceptSpec *tmp = pmConceptSpecAlloc(blank, dummyConceptParser, - dummyConceptFormatter, true); + dummyConceptFormatter, dummyConceptCopier, true); ok(tmp != NULL, "pmConceptSpecAlloc() returned non-NULL"); skip_start(tmp == NULL, 4, "Skipping tests because pmConceptSpecAlloc() returned NULL"); @@ -98,4 +116,5 @@ ok(tmp->parse == dummyConceptParser, "pmConceptSpecAlloc() set the ->parse member correctly"); ok(tmp->format == dummyConceptFormatter, "pmConceptSpecAlloc() set the ->format member correctly"); + ok(tmp->copy == dummyConceptCopier, "pmConceptSpecAlloc() set the ->copy member correctly"); ok(tmp->required == true, "pmConceptSpecAlloc() set the ->required member correctly"); skip_end(); @@ -109,5 +128,5 @@ { psMemId id = psMemGetId(); - pmConceptSpec *tmp = pmConceptSpecAlloc(NULL, NULL, NULL, false); + pmConceptSpec *tmp = pmConceptSpecAlloc(NULL, NULL, NULL, NULL, false); ok(tmp != NULL, "pmConceptSpecAlloc() returned non-NULL with NULL inputs"); psFree(tmp); Index: branches/pap/psModules/test/objects/tap_pmGrowthCurve.c =================================================================== --- branches/pap/psModules/test/objects/tap_pmGrowthCurve.c (revision 23948) +++ branches/pap/psModules/test/objects/tap_pmGrowthCurve.c (revision 25027) @@ -129,5 +129,5 @@ source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1; - source->mode = PM_SOURCE_MODE_SUBTRACTED; + source->mode = PM_SOURCE_MODE_PSFSTAR; source->modelPSF->radiusFit = 15.0; @@ -232,5 +232,5 @@ source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1; - source->mode = PM_SOURCE_MODE_SUBTRACTED; + source->mode = PM_SOURCE_MODE_PSFSTAR; source->modelPSF->radiusFit = 15.0; @@ -334,5 +334,5 @@ source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1; - source->mode = PM_SOURCE_MODE_SUBTRACTED; + source->mode = PM_SOURCE_MODE_PSFSTAR; source->modelPSF->radiusFit = 15.0; Index: branches/pap/psModules/test/objects/tap_pmSource.c =================================================================== --- branches/pap/psModules/test/objects/tap_pmSource.c (revision 23948) +++ branches/pap/psModules/test/objects/tap_pmSource.c (revision 25027) @@ -432,5 +432,5 @@ psArray *sources = psArrayAlloc(NUM_SOURCES); psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); @@ -452,5 +452,5 @@ psArray *sources = psArrayAlloc(NUM_SOURCES); psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); @@ -490,5 +490,5 @@ } psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); Index: branches/pap/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c =================================================================== --- branches/pap/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c (revision 23948) +++ branches/pap/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c (revision 25027) @@ -48,5 +48,5 @@ psMetadata *tableHeader = psMetadataAlloc(); psString extname = psStringCopy("ext"); - bool rc = pmSourcesWrite_PS1_DEV_1(NULL, sources, imageHeader, tableHeader, extname, NULL); + bool rc = pmSourcesWrite_PS1_DEV_1(NULL, sources, imageHeader, tableHeader, extname); ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL psFits input parameter"); psFree(fitsFile); @@ -72,5 +72,5 @@ psMetadata *tableHeader = psMetadataAlloc(); psString extname = psStringCopy("ext"); - bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, NULL, imageHeader, tableHeader, extname, NULL); + bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, NULL, imageHeader, tableHeader, extname); ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL pmSource input parameter"); psFree(fitsFile); @@ -96,5 +96,5 @@ psMetadata *tableHeader = psMetadataAlloc(); psString extname = psStringCopy("ext"); - bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, NULL, NULL); + bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, NULL); ok(rc == false, "pmSourcesWrite_PS1_DEV_1() returned FALSE with NULL extname input parameter"); psFree(fitsFile); @@ -216,5 +216,5 @@ psMetadata *tableHeader = psMetadataAlloc(); psString extname = psStringCopy("ext"); - bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, extname, NULL); + bool rc = pmSourcesWrite_PS1_DEV_1(fitsFile, sources, imageHeader, tableHeader, extname); ok(rc == true, "pmSourcesWrite_PS1_DEV_1() returned TRUE with acceptable input parameters"); psFree(fitsFile);