Index: /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.c (revision 25407) @@ -289,4 +289,5 @@ // test the CDELTi varient if (pcKeys) { + wcs->wcsCDkeys = 0; wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1"); wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2"); @@ -334,4 +335,6 @@ // test the CDi_j varient if (cdKeys) { + wcs->wcsCDkeys = 1; + wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1 wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2 @@ -375,42 +378,59 @@ // XXX make it optional to write out CDi_j terms, or other versions // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity - double cdelt1 = wcs->cdelt1; - double cdelt2 = wcs->cdelt2; - psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); - psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); - - // test the PC00i00j varient: - psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1 - psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2 - psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1 - psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2 - - // Elixir-style polynomial terms - // XXX currently, Elixir/DVO cannot accept mixed orders - // XXX need to respect the masks - // XXX is wcs->cdelt1,2 always consistent? - int fitOrder = wcs->trans->x->nX; - if (fitOrder > 1) { + if (!(wcs->wcsCDkeys)) { + + double cdelt1 = wcs->cdelt1; + double cdelt2 = wcs->cdelt2; + psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); + psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); + + // test the PC00i00j varient: + psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1 + psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2 + psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1 + psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2 + + // Elixir-style polynomial terms + // XXX currently, Elixir/DVO cannot accept mixed orders + // XXX need to respect the masks + // XXX is wcs->cdelt1,2 always consistent? + int fitOrder = wcs->trans->x->nX; + if (fitOrder > 1) { for (int i = 0; i <= fitOrder; i++) { - for (int j = 0; j <= fitOrder; j++) { - if (i + j < 2) - continue; - if (i + j > fitOrder) - continue; - sprintf (name, "PCA1X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - sprintf (name, "PCA2X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - } + for (int j = 0; j <= fitOrder; j++) { + if (i + j < 2) + continue; + if (i + j > fitOrder) + continue; + sprintf (name, "PCA1X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + sprintf (name, "PCA2X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + } } psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder); - } - - // remove any existing 'CDi_j style' wcs keywords - if (psMetadataLookup(header, "CD1_1")) { + } + + // remove any existing 'CDi_j style' wcs keywords + if (psMetadataLookup(header, "CD1_1")) { psMetadataRemoveKey(header, "CD1_1"); psMetadataRemoveKey(header, "CD1_2"); psMetadataRemoveKey(header, "CD2_1"); psMetadataRemoveKey(header, "CD2_2"); + } + } + if (wcs->wcsCDkeys) { + + psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]); + psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]); + psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]); + psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]); + + if (psMetadataLookup(header, "PC001001")) { + psMetadataRemoveKey(header, "PC001001"); + psMetadataRemoveKey(header, "PC001002"); + psMetadataRemoveKey(header, "PC002001"); + psMetadataRemoveKey(header, "PC002002"); + } } @@ -535,4 +555,6 @@ fpa->toSky->R -= 2.0*M_PI; + fpa->wcsCDkeys = wcs->wcsCDkeys; + psTrace ("psastro", 5, "toFPA: %f %f (%f,%f),(%f,%f)\n", chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], @@ -682,4 +704,5 @@ wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]); + wcs->wcsCDkeys = fpa->wcsCDkeys; psFree (toTPA); @@ -922,4 +945,5 @@ wcs->trans = psPlaneTransformAlloc (nXorder, nYorder); wcs->toSky = NULL; + wcs->wcsCDkeys = 0; memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE); Index: /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryWCS.h (revision 25407) @@ -23,4 +23,5 @@ double crpix1, crpix2; double cdelt1, cdelt2; + bool wcsCDkeys; psProjection *toSky; psPlaneTransform *trans; Index: /branches/eam_branches/20090715/psModules/src/camera/pmFPA.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmFPA.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmFPA.h (revision 25407) @@ -48,4 +48,5 @@ psPlaneTransform *toTPA; ///< Transformation from focal plane to tangent plane, or NULL psProjection *toSky; ///< Projection from tangent plane to sky, or NULL + bool wcsCDkeys; // Information psMetadata *concepts; ///< FPA-level concepts Index: /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.c (revision 25407) @@ -111,6 +111,6 @@ // 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; + psWarning("CELL.SATURATION is not set --- completely masking cell.\n"); + saturation = NAN; } float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level @@ -118,10 +118,10 @@ // 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; + 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; + // 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 @@ -140,6 +140,6 @@ // completely mask if SATURATION or BAD are invalid if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) { - psImageInit(mask, badMask); - return true; + psImageInit(mask, badMask); + return true; } @@ -230,5 +230,5 @@ // return false; psWarning("CELL.GAIN is not set --- setting variance to NAN\n"); - gain = NAN; + gain = NAN; } float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise @@ -237,5 +237,5 @@ // return false; psWarning("CELL.READNOISE is not set --- setting variance to NAN\n"); - readnoise = NAN; + readnoise = NAN; } // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN) @@ -248,10 +248,10 @@ if (isnan(gain) || isnan(readnoise)) { if (!readout->variance) { - // generate the image if needed + // 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 + // XXX need to set the mask, if defined psImageInit(readout->variance, NAN); - return true; + return true; } @@ -262,5 +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 + // 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", @@ -276,9 +276,9 @@ // 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); + 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)); + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); } @@ -362,6 +362,6 @@ -bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal, - psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) +bool pmReadoutVarianceRenormalise(const pmReadout *readout, psImageMaskType maskVal, + int sample, float minValid, float maxValid) { PM_ASSERT_READOUT_NON_NULL(readout, false); @@ -371,263 +371,80 @@ psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts - if (!psMemIncrRefCounter(rng)) { - rng = psRandomAlloc(PS_RANDOM_TAUS); - } - - psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean - if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image"); - psFree(varianceStats); - psFree(rng); - return false; - } - float meanVariance = varianceStats->robustMedian; // Mean variance - psFree(varianceStats); - - psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean - if (!psImageBackground(imageStats, NULL, image, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure stdev of image"); - psFree(imageStats); - psFree(rng); - return false; - } - float stdevImage = imageStats->robustStdev; // Standard deviation of image - psFree(imageStats); - psFree(rng); - - float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be - psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction); - psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32)); - - return true; -} - - -bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width, - psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng) -{ - PM_ASSERT_READOUT_NON_NULL(readout, false); - PM_ASSERT_READOUT_IMAGE(readout, false); - PM_ASSERT_READOUT_VARIANCE(readout, false); - - if (!psMemIncrRefCounter(rng)) { - rng = psRandomAlloc(PS_RANDOM_TAUS); - } - - psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images - - // Measure background - psStats *bgStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);// Statistics for background - if (!psImageBackground(bgStats, NULL, image, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image"); - psFree(bgStats); - psFree(rng); - return false; - } - float bgMean = bgStats->robustMedian; // Background level - float bgNoise = bgStats->robustStdev; // Background standard deviation - psFree(bgStats); - psTrace("psModules.camera", 5, "Background is %f +/- %f\n", bgMean, bgNoise); - - - // Construct kernels for flux measurement - // We use N(0,width) and N(0,width/sqrt(2)) kernels, following psphotSignificanceImage. - float sigFactor = 4.0 * M_PI * PS_SQR(width); // Factor for conversion from im/wt ratio to significance - int size = RENORM_NUM_SIGMA, fullSize = 2 * size + 1; // Half-size and full size of Gaussian - psVector *gauss = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian for weighting - psVector *gauss2 = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian squared - for (int i = 0, x = -size; i < fullSize; i++, x++) { - gauss->data.F32[i] = expf(-0.5 * PS_SQR(x) / PS_SQR(width)); - gauss2->data.F32[i] = expf(-PS_SQR(x) / PS_SQR(width)); - } - - // Size of image - int numCols = image->numCols, numRows = image->numRows; // Size of images - int xSize = numCols - fullSize, ySize = numRows - fullSize; // Size of consideration - int xOffset = size, yOffset = size; // Offset to region of consideration - - // Measure fluxes - float peakFlux = RENORM_PEAK * bgNoise; // Peak flux for fake sources - psVector *noise = psVectorAlloc(num, PS_TYPE_F32); // Measurements of the noise - psVector *source = psVectorAlloc(num, PS_TYPE_F32); // Measurements of fake sources - psVector *guess = psVectorAlloc(num, PS_TYPE_F32); // Guess at significance - psVector *photMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for fluxes - for (int i = 0; i < num; i++) { - // Coordinates of interest - int xPix = psRandomUniform(rng) * xSize + xOffset + 0.5; - int yPix = psRandomUniform(rng) * ySize + yOffset + 0.5; - psAssert(xPix - size >= 0 && xPix + size < numCols && - yPix - size >= 0 && yPix + size < numRows, - "Bad pixel position: %d,%d", xPix, yPix); - - // Weighted aperture photometry - // This has the same effect as smoothing the image by the window function - float sumNoise = 0.0; // Sum for noise measurement - float sumSource = 0.0; // Sum for source measurement - float sumVariance = 0.0; // Sum for variance measurement - float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels - for (int v = 0, y = yPix - size; v < fullSize; v++, y++) { - float xSumNoise = 0.0; // Sum for noise measurement in x - float xSumSource = 0.0; // Sum for source measurement in x - float xSumVariance = 0.0; // Sum for variance measurement in x - float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x - float yGauss = gauss->data.F32[v]; // Value of Gaussian in y - for (int u = 0, x = xPix - size; u < fullSize; u++, x++) { - if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { + int numCols = image->numCols, numRows = image->numRows; // Size of image + int numPix = numCols * numRows; // Number of pixels + int num = PS_MAX(sample, numPix); // Number we care about + psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32); // Signal-to-noise values + + if (num >= numPix) { + // We have an image smaller than Nsubset, so just loop over the image pixels + int index = 0; // Index for vector + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || + !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) { continue; } - float value = image->data.F32[y][x] - bgMean; // Value of image - float xGauss = gauss->data.F32[u]; // Value of Gaussian in x - float xGauss2 = gauss2->data.F32[u]; // Value of Gaussian^2 in x - xSumNoise += value * xGauss; - xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss; - xSumVariance += variance->data.F32[y][x] * xGauss2; - xSumGauss += xGauss; - xSumGauss2 += xGauss2; - } - float yGauss2 = gauss2->data.F32[v]; // Value of Gaussian^2 in y - sumNoise += xSumNoise * yGauss; - sumSource += xSumSource * yGauss; - sumVariance += xSumVariance * yGauss2; - sumGauss += xSumGauss * yGauss; - sumGauss2 += xSumGauss2 * yGauss2; - } - - photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) && - isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ? - 0 : 0xFF); - - float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise - float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source - float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel - - noise->data.F32[i] = smoothImageNoise; - source->data.F32[i] = smoothImageSource; - guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0; - psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n", - i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance); - } - psFree(gauss); - psFree(gauss2); - psFree(rng); - - // Standard deviation of fluxes gives us the real significance - psStats *stdevStats = psStatsAlloc(stdevStat); // Statistics - if (!psVectorStats(stdevStats, noise, NULL, photMask, 0xFF)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure standard deviation of fluxes"); - psFree(stdevStats); - psFree(noise); - psFree(source); - psFree(guess); - psFree(photMask); + + signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]); + index++; + } + } + signoise->n = index; + } else { + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + int index = 0; // Index for vector + for (long i = 0; i < num; i++) { + // Pixel coordinates + int pixel = numPix * psRandomUniform(rng); + int x = pixel % numCols; + int y = pixel / numCols; + + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || + !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) { + continue; + } + + signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]); + index++; + } + signoise->n = index; + psFree(rng); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics + + if (!psVectorStats(stats, signoise, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics on S/N image"); + psFree(signoise); return false; } - float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes - psFree(stdevStats); - psFree(noise); - psTrace("psModules.camera", 5, "Standard deviation of fluxes is %f\n", stdev); - - // Ratio of measured significance to guessed significance - psVector *ratio = psVectorAlloc(num, PS_TYPE_F32); // Ratio of measured to guess - for (int i = 0; i < num; i++) { - float measuredSig = PS_SQR(source->data.F32[i] / stdev); // Measured significance - ratio->data.F32[i] = measuredSig / guess->data.F32[i]; - if (guess->data.F32[i] <= 0.0 || source->data.F32[i] <= 0.0 || !isfinite(ratio->data.F32[i])) { - photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF; - } - psTrace("psModules.camera", 9, "Ratio %d: %f, %f, %f\n", - i, guess->data.F32[i], measuredSig, ratio->data.F32[i]); - } - psFree(source); - psFree(guess); - - psStats *meanStats = psStatsAlloc(meanStat | stdevStat); // Statistics - if (!psVectorStats(meanStats, ratio, NULL, photMask, 0xFF)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure mean ratio"); - psFree(meanStats); - psFree(ratio); - psFree(photMask); - return false; - } - float meanRatio = psStatsGetValue(meanStats, meanStat); // Mean ratio - psTrace("psModules.camera", 5, "Mean significance ratio is %f +/- %f\n", - meanRatio, psStatsGetValue(meanStats, stdevStat)); - psFree(meanStats); - psFree(ratio); - psFree(photMask); - - psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio); - psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32)); + psFree(signoise); + + float covar = sqrtf(psImageCovarianceFactor(readout->covariance)); // Covariance factor + float correction = stats->robustStdev / covar; // Correction factor + psFree(stats); + psLogMsg("psModules.camera", PS_LOG_DETAIL, "Variance renormalisation factor is %f", correction); + + // Check valid range of correction factor + if ((isfinite(minValid) && correction < minValid) || (isfinite(maxValid) && correction > maxValid)) { + psWarning("Variance renormalisation is outside valid range: %f vs %f:%f --- no correction made", + correction, minValid, maxValid); + return true; + } + + psBinaryOp(variance, variance, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32)); + + pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout + if (hdu) { + psString history = NULL; + psStringAppend(&history, "Rescaled variance by %6.4f (stdev by %6.4f)", + PS_SQR(correction), correction); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + psFree(history); + } return true; } - -bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat, - psStatsOptions stdevStat, int width, psRandom *rng) -{ - PM_ASSERT_READOUT_NON_NULL(readout, false); - PM_ASSERT_READOUT_IMAGE(readout, false); - PM_ASSERT_READOUT_VARIANCE(readout, false); - PS_ASSERT_INT_POSITIVE(width, false); - - if (!psMemIncrRefCounter(rng)) { - rng = psRandomAlloc(PS_RANDOM_TAUS); - } - - psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images - int numCols = image->numCols, numRows = image->numRows; // Size of images - int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions - float xSize = numCols / (float)xNum, ySize = numRows / (float)yNum; // Size of renormalisation regions - - psStats *meanStats = psStatsAlloc(meanStat), *stdevStats = psStatsAlloc(stdevStat); // Statistics - psVector *buffer = NULL; - - for (int j = 0; j < yNum; j++) { - // Bounds in y - int yMin = j * ySize; - int yMax = (j + 1) * ySize; - for (int i = 0; i < xNum; i++) { - // Bounds in x - int xMin = i * xSize; - int xMax = (i + 1) * xSize; - - psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest - psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels - psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels - psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels - - if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) || - !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) { - // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it. - psString regionStr = psRegionToString(region); // String with region - psWarning("Unable to measure statistics over %s", regionStr); - psFree(regionStr); - psErrorClear(); - psImageInit(subVariance, NAN); - if (subMask) { - psImageInit(subMask, maskVal); - } - } else { - double meanVar = psStatsGetValue(meanStats, meanStat); // Mean of variance map - double stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of image - psTrace("psModules.camera", 3, - "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n", - xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar); - psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32)); - } - - psFree(subImage); - psFree(subVariance); - psFree(subMask); - } - } - psFree(meanStats); - psFree(stdevStats); - psFree(rng); - psFree(buffer); - - return true; -} Index: /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmFPAMaskWeight.h (revision 25407) @@ -15,24 +15,4 @@ /// @addtogroup Camera Camera Layout /// @{ - -#if 0 -/// Pixel mask values -typedef enum { - PM_MASK_CLEAR = 0x00, ///< The pixel is not masked - PM_MASK_BLANK = 0x01, ///< The pixel is blank or has no (valid) data - PM_MASK_FLAT = 0x02, ///< The pixel is non-positive in the flat-field - PM_MASK_DETECTOR = 0x02, ///< The detector pixel is bad (e.g., bad column, charge trap) - PM_MASK_SAT = 0x04, ///< The pixel is saturated in the image of interest - PM_MASK_BAD = 0x04, ///< The pixel is low in the image of interest - PM_MASK_RANGE = 0x04, ///< The pixel is out of range in the image of interest - PM_MASK_CR = 0x08, ///< The pixel is probably a CR - PM_MASK_SPARE1 = 0x10, ///< Spare mask value - PM_MASK_SPARE2 = 0x20, ///< Spare mask value - PM_MASK_SUSPECT = 0x40, ///< The pixel is suspected of being bad, but may not be - PM_MASK_MARK = 0x80, ///< The pixel is marked as temporarily ignored -} pmMaskValue; -#define PM_MASK_MARK 0x80 ///< The pixel is marked as temporarily ignored -#define PM_MASK_SAT 0x04 ///< The pixel is saturated in the image of interest -#endif /// Set a temporary readout mask using CELL.SATURATION and CELL.BAD @@ -54,5 +34,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 + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -73,5 +53,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 + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -83,5 +63,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 + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -93,5 +73,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 + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -100,36 +80,10 @@ /// /// The variance map is adjusted so that the mean matches the actual pixel variance in the image -bool pmReadoutVarianceRenormPixels( - const pmReadout *readout, ///< Readout to normalise - psImageMaskType maskVal, ///< Value to mask - psStatsOptions meanStat, ///< Statistic to measure the mean (of the variance map) - psStatsOptions stdevStat, ///< Statistic to measure the stdev (of the image) - psRandom *rng ///< Random number generator - ); - -/// Renormalise the variance map to match the actual photometry variance -/// -/// The variance map is adjusted so that the actual significance of fake sources matches the -/// guestimated significance -bool pmReadoutVarianceRenormPhot( +bool pmReadoutVarianceRenormalise( const pmReadout *readout, ///< Readout to normalise psImageMaskType maskVal, ///< Value to mask - int num, ///< Number of instances to measure over the image - float width, ///< Photometry width - psStatsOptions meanStat, ///< Statistic to measure the mean - psStatsOptions stdevStat, ///< Statistic to measure the stdev - psRandom *rng ///< Random number generator - ); - -/// Renormalise the variance map to match the actual variance -/// -/// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for -/// that patch corresponds. -bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise - psImageMaskType maskVal, // Value to mask - psStatsOptions meanStat, // Statistic to measure the mean (of the variance map) - psStatsOptions stdevStat, // Statistic to measure the stdev (of the image) - int width, // Width of patch (pixels) - psRandom *rng // Random number generator (for sub-sampling images) + int sample, ///< Sample size + float minValid, ///< Minimum valid renormalisation, or NAN + float maxValid ///< Maximum valid renormalisation, or NAN ); Index: /branches/eam_branches/20090715/psModules/src/camera/pmFPAMosaic.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmFPAMosaic.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmFPAMosaic.c (revision 25407) @@ -740,5 +740,5 @@ static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned psImage **mosaicMask, // The mosaic mask, to be returned - psImage **mosaicVariances, // The mosaic variances, to be returned + psImage **mosaicVariance, // The mosaic variance, to be returned int *xBinChip, int *yBinChip, // The binning in x and y, to be returned const pmChip *chip, // Chip to mosaic @@ -749,5 +749,5 @@ assert(mosaicImage); assert(mosaicMask); - assert(mosaicVariances); + assert(mosaicVariance); assert(xBinChip); assert(yBinChip); @@ -826,5 +826,5 @@ if (allGood) { *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); - *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); + *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE); *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank); } @@ -847,5 +847,5 @@ static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned psImage **mosaicMask, // The mosaic mask, to be returned - psImage **mosaicVariances, // The mosaic variances, to be returned + psImage **mosaicVariance, // The mosaic variance, to be returned int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned const pmFPA *fpa, // FPA to mosaic @@ -857,5 +857,5 @@ assert(mosaicImage); assert(mosaicMask); - assert(mosaicVariances); + assert(mosaicVariance); assert(xBinFPA); assert(yBinFPA); @@ -960,5 +960,5 @@ if (allGood) { *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); - *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); + *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE); *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank); } @@ -1025,5 +1025,5 @@ psImage *mosaicImage = NULL; // The mosaic image psImage *mosaicMask = NULL; // The mosaic mask - psImage *mosaicVariances = NULL; // The mosaic variances + psImage *mosaicVariance = NULL; // The mosaic variances // Find the HDU @@ -1052,6 +1052,6 @@ } if (hdu->variances) { - mosaicVariances = psImageSubset(hdu->variances->data[0], bounds); - if (!mosaicVariances) { + mosaicVariance = psImageSubset(hdu->variances->data[0], bounds); + if (!mosaicVariance) { psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n"); return false; @@ -1061,5 +1061,5 @@ // Case 2 --- we need to mosaic by cut and paste psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n"); - if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) { + if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariance, &xBin, &yBin, source, targetCell, blank)) { psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n"); return false; @@ -1069,4 +1069,5 @@ } psTrace("psModules.camera", 1, "xBin,yBin: %d,%d\n", xBin, yBin); + // Set the concepts for the target cell @@ -1090,9 +1091,31 @@ target->parent->concepts = psMetadataCopy(target->parent->concepts, source->parent->concepts); // FPA lvl + // Average the covariances + psList *covariances = psListAlloc(NULL); // Input covariance matrices + for (int i = 0; i < source->cells->n; i++) { + pmCell *cell = source->cells->data[i]; // Cell of interest + if (!cell || !cell->data_exists) { + continue; + } + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + if (!ro || !ro->covariance) { + continue; + } + psListAdd(covariances, PS_LIST_TAIL, ro->covariance); + } + psKernel *mosaicCovariance = NULL; // Covariance for mosaic + if (psListLength(covariances) > 0) { + psArray *covarArray = psListToArray(covariances); // Array with covariances + mosaicCovariance = psImageCovarianceAverage(covarArray); + psFree(covarArray); + } + psFree(covariances); + // Now make a new readout to go in the target cell pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout newReadout->image = mosaicImage; newReadout->mask = mosaicMask; - newReadout->variance = mosaicVariances; + newReadout->variance = mosaicVariance; + newReadout->covariance = mosaicCovariance; psFree(newReadout); // Drop reference @@ -1334,4 +1357,32 @@ target->concepts = psMetadataCopy(target->concepts, source->concepts); + // Average the covariances + psList *covariances = psListAlloc(NULL); // Input covariance matrices + for (int i = 0; i < covariances->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + if (!chip || !chip->data_exists) { + continue; + } + psArray *cells = chip->cells; // Cells in chip + for (long j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[i]; // Cell of interest + if (!cell || !cell->data_exists) { + continue; + } + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + if (!ro || !ro->covariance) { + continue; + } + psListAdd(covariances, PS_LIST_TAIL, ro->covariance); + } + } + psKernel *mosaicCovariances = NULL; // Covariance for mosaic + if (psListLength(covariances) > 0) { + psArray *covarArray = psListToArray(covariances); // Array with covariances + mosaicCovariances = psImageCovarianceAverage(covarArray); + psFree(covarArray); + } + psFree(covariances); + // Now make a new readout to go in the new cell pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout @@ -1339,4 +1390,5 @@ newReadout->mask = mosaicMask; newReadout->variance = mosaicVariances; + newReadout->covariance = mosaicCovariances; psFree(newReadout); // Drop reference Index: /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.c (revision 25407) @@ -28,5 +28,4 @@ #define MODEL_TYPE "PS_MODEL_RGAUSS" // Type of model to use #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 @@ -50,13 +49,23 @@ } -bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources, - const psVector *xOffset, const psVector *yOffset, const pmPSF *psf, - float minFlux, int radius, bool circularise, bool normalisePeak) + +bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows, + const psVector *x, const psVector *y, const psVector *mag, + const psVector *xOffset, const psVector *yOffset, + const pmPSF *psf, float minFlux, int radius, + bool circularise, bool normalisePeak) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_INT_LARGER_THAN(numCols, 0, false); PS_ASSERT_INT_LARGER_THAN(numRows, 0, false); - PS_ASSERT_ARRAY_NON_NULL(sources, false); - + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false); + PS_ASSERT_VECTOR_NON_NULL(mag, false); + PS_ASSERT_VECTOR_TYPE(mag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(mag, x, false); + long numSources = x->n; // Number of sources if (xOffset || yOffset) { PS_ASSERT_VECTOR_NON_NULL(xOffset, false); @@ -64,47 +73,24 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false); PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false); - if (xOffset->n != sources->n) { + PS_ASSERT_VECTOR_TYPE(yOffset, PS_TYPE_S32, false); + if (xOffset->n != numSources) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of offset vectors (%ld) and sources (%ld) doesn't match", - xOffset->n, sources->n); + xOffset->n, numSources); return false; } } PS_ASSERT_PTR_NON_NULL(psf, false); - if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius."); - return false; - } readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32); psImageInit(readout->image, 0); - int numSources = sources->n; // Number of stars - for (int 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; - } - float x, y; // Coordinates of source - if (source->modelPSF) { - x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - } else { - x = source->peak->xf; - y = source->peak->yf; - } - - float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source + for (long i = 0; i < numSources; i++) { + float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source + float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source if (normalisePeak) { // Normalise flux - pmModel *normModel = pmModelFromPSFforXY(psf, x, y, 1.0); // Model for normalisation + pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation if (!normModel || (normModel->flags & MODEL_MASK)) { psFree(normModel); @@ -121,5 +107,5 @@ } - pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, flux); + pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux); if (!fakeModel || (fakeModel->flags & MODEL_MASK)) { psFree(fakeModel); @@ -137,11 +123,16 @@ pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate - fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); - float fakeRadius = radius > 0 ? radius : - PS_MAX(1.0, fakeModel->modelRadius(fakeModel->params, minFlux)); // Radius of fake source + fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); + float fakeRadius = 1.0; // Radius of fake source + if (isfinite(minFlux)) { + fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux)); + } + if (radius > 0) { + fakeRadius = PS_MAX(fakeRadius, radius); + } if (xOffset) { - if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i], - y + yOffset->data.S32[i], fakeRadius)) { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i], + ySrc + yOffset->data.S32[i], fakeRadius)) { psErrorClear(); continue; @@ -153,5 +144,5 @@ } } else { - if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) { psErrorClear(); continue; @@ -168,2 +159,54 @@ return true; } + + +bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources, + pmSourceMode sourceMask, const psVector *xOffset, const psVector *yOffset, + const pmPSF *psf, float minFlux, int radius, + bool circularise, bool normalisePeak) +{ + PS_ASSERT_ARRAY_NON_NULL(sources, false); + + int numSources = sources->n; // Number of stars + psVector *x = psVectorAllocEmpty(numSources, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty(numSources, PS_TYPE_F32); + psVector *mag = psVectorAllocEmpty(numSources, PS_TYPE_F32); + + int numGood = 0; // Number of good sources + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + if (source->mode & sourceMask) { + continue; + } + if (!isfinite(source->psfMag)) { + continue; + } + float xSrc, ySrc; // Coordinates of source + if (source->modelPSF) { + xSrc = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + ySrc = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + xSrc = source->peak->xf; + ySrc = source->peak->yf; + } + + x->data.F32[numGood] = xSrc; + y->data.F32[numGood] = ySrc; + mag->data.F32[numGood] = source->psfMag; + numGood++; + } + x->n = numGood; + y->n = numGood; + mag->n = numGood; + + bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf, + minFlux, radius, circularise, normalisePeak); + psFree(x); + psFree(y); + psFree(mag); + + return status; +} Index: /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.h (revision 25407) @@ -11,9 +11,25 @@ #include #include +#include + +/// Generate a fake readout from vectors +bool pmReadoutFakeFromVectors(pmReadout *readout, ///< Output readout + int numCols, int numRows, ///< Dimension of image + const psVector *x, const psVector *y, ///< Source coordinates + const psVector *mag, ///< Source magnitudes + const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL + const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL + const pmPSF *psf, ///< PSF for sources + float minFlux, ///< Minimum flux to bother about; for setting source radius + int radius, ///< Fixed radius for sources + bool circularise, ///< Circularise PSF model? + bool normalisePeak ///< Normalise the peak value? + ); /// Generate a fake readout from an array of sources -bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL +bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout int numCols, int numRows, ///< Dimension of image const psArray *sources, ///< Array of pmSource + pmSourceMode sourceMask, ///< Mask for sources const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL Index: /branches/eam_branches/20090715/psModules/src/concepts/pmConcepts.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/concepts/pmConcepts.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/concepts/pmConcepts.c (revision 25407) @@ -377,5 +377,4 @@ conceptRegisterS32("CELL.XWINDOW", "Start of cell window (pixels)",p_pmConceptParse_Positions,p_pmConceptFormat_Positions, NULL, true, PM_FPA_LEVEL_CELL); conceptRegisterS32("CELL.YWINDOW", "Start of cell window (pixels)",p_pmConceptParse_Positions,p_pmConceptFormat_Positions, NULL, true, PM_FPA_LEVEL_CELL); - conceptRegisterF32("CELL.VARFACTOR", "Variance factor for conversion from large to small scales", NULL, NULL, NULL, true, PM_FPA_LEVEL_CELL); // CELL.TRIMSEC @@ -485,27 +484,27 @@ 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; - } + // 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/eam_branches/20090715/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/concepts/pmConceptsStandard.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/concepts/pmConceptsStandard.c (revision 25407) @@ -429,8 +429,10 @@ } - 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; + if (strcmp(concept->name, "FPA.RA") == 0 || strcmp(concept->name, "FPA.DEC") == 0) { + 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 { @@ -450,9 +452,8 @@ 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); + psStringAppend(&coordString, "%s%02d:%02d:%06.3f", + negative ? "-" : (strcmp(concept->name, "FPA.DEC") == 0 ? "+" : ""), + big, medium, small); coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString); psFree(coordString); Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c (revision 25407) @@ -124,5 +124,5 @@ pmResiduals *resid = psf->residuals;// PSF residuals psf->residuals = NULL; - if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf, + if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf, NAN, radius, true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout."); @@ -298,5 +298,5 @@ } - // measure the source moments: tophat windowing, no pixel S/N cutoff + // 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 Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmStack.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmStack.c (revision 25407) @@ -30,5 +30,5 @@ #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time -#define ADD_VARIANCE // Allow additional variance (besides variance factor)? +//#define ADD_VARIANCE // Allow additional variance (besides variance factor)? #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.c (revision 25407) @@ -733,5 +733,5 @@ int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, - const psVector *deviations, psImage *subMask, float sigmaRej, int footprint) + const psVector *deviations, psImage *subMask, float sigmaRej) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); @@ -821,4 +821,5 @@ ds9num++; + int footprint = stamps->footprint; // Half-size of stamp region of interest int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps @@ -956,10 +957,10 @@ PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return - psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector + psArray *images = psArrayAlloc(kernels->solution1->n - 1); // Images of each kernel to return + psVector *fakeSolution = psVectorAlloc(kernels->solution1->n, PS_TYPE_F64); // Fake solution vector psVectorInit(fakeSolution, 0.0); - for (int i = 0; i < solution->n - 1; i++) { - fakeSolution->data.F64[i] = solution->data.F64[i]; + for (int i = 0; i < kernels->solution1->n - 1; i++) { + fakeSolution->data.F64[i] = kernels->solution1->data.F64[i]; images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); fakeSolution->data.F64[i] = 0.0; Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtraction.h (revision 25407) @@ -68,6 +68,5 @@ const psVector *deviations, ///< Deviations for each stamp psImage *subMask, ///< Subtraction mask - float sigmaRej, ///< Number of RMS deviations above zero at which to reject - int footprint ///< Half-size of stamp + float sigmaRej ///< Number of RMS deviations above zero at which to reject ); @@ -94,8 +93,9 @@ /// Generate images of the convolution kernel elements -psArray *pmSubtractionKernelSolutions(const psVector *solution, ///< Solution vector - const pmSubtractionKernels *kernels, ///< Kernel parameters - float x, float y ///< Normalised position [-1,1] for images +psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y, ///< Normalised position [-1,1] for images + bool wantDual ///< Calculate for the dual kernel? ); + /// Execute a thread job to convolve a patch of the image Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25407) @@ -17,10 +17,10 @@ -bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region, +bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, + pmSubtractionKernels *kernels, psRegion *region, int numCols, int numRows) { - if (analysis) { - PS_ASSERT_METADATA_NON_NULL(analysis, false); - } + PS_ASSERT_METADATA_NON_NULL(analysis, false); + PS_ASSERT_METADATA_NON_NULL(header, false); PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); @@ -38,5 +38,11 @@ PS_DATA_REGION | PS_META_DUPLICATE_OK, "Region over which subtraction was performed", subRegion); + + psString string = psRegionToString(*subRegion); psFree(subRegion); + + psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", string); + psFree(string); } @@ -45,5 +51,7 @@ PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, - PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->mode); + PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode); + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode); // Realisations of kernel @@ -113,11 +121,11 @@ { psMetadata *header = psMetadataAlloc(); // Header - for (int i = 0; i < solution->n; i++) { + for (int i = 0; i < kernels->solution1->n; i++) { psString name = NULL; // Header keyword psStringAppend(&name, "SOLN%04d", i); - psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]); + psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution1->data.F64[i]); psFree(name); } - psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0); + psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false); psFits *kernelFile = psFitsOpen("kernels.fits", "w"); (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); @@ -128,42 +136,4 @@ #endif - - // Set the variance factors - float vf1 = 1.0, vf2 = 1.0; // Variance factors for each image - switch (kernels->mode) { - case PM_SUBTRACTION_MODE_1: - vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); - break; - case PM_SUBTRACTION_MODE_2: - vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); - break; - case PM_SUBTRACTION_MODE_DUAL: - vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false); - vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true); - break; - default: - psAbort("Invalid subtraction mode: %x", kernels->mode); - } - - // Weight by the area - if (region) { - float norm = (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) / (numCols * numRows); - vf1 *= norm; - vf2 *= norm; - } - - // Update the variance factor -#define UPDATE_VARFACTOR(VF, ANALYSIS) { \ - psMetadataItem *vfItem = psMetadataLookup(analysis, ANALYSIS); \ - if (vfItem) { \ - psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said."); \ - vfItem->data.F32 += VF; \ - } else { \ - psMetadataAddF32(analysis, PS_LIST_TAIL, ANALYSIS, 0, "Variance factor weighted by the area", VF); \ - } \ -} - - UPDATE_VARFACTOR(vf1, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1); - UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2); // Kernel shape @@ -201,10 +171,22 @@ psFree(image); - psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous - if (item) { - item->data.F32 = PS_MAX(item->data.F32, max); - } else { - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, - "Maximum deconvolution fraction", max); + { + psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + max = item->data.F32 = PS_MAX(item->data.F32, max); + } else { + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, + "Maximum deconvolution fraction", max); + } + } + + { + psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + item->data.F32 = max; + } else { + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, + "Maximum deconvolution fraction", max); + } } } @@ -254,4 +236,17 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY, PS_META_DUPLICATE_OK, "Moment in yy", m02); + + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, + PS_META_DUPLICATE_OK, "Normalisation", m00); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX, + PS_META_DUPLICATE_OK, "Moment in x", m10); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY, + PS_META_DUPLICATE_OK, "Moment in y", m01); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX, + PS_META_DUPLICATE_OK, "Moment in xx", m20); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY, + PS_META_DUPLICATE_OK, "Moment in xy", m11); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY, + PS_META_DUPLICATE_OK, "Moment in yy", m02); } @@ -263,4 +258,6 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF, + PS_META_DUPLICATE_OK, "Background difference", bg); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF, PS_META_DUPLICATE_OK, "Background difference", bg); psFree(polyValues); @@ -275,4 +272,11 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", kernels->rms); + + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps", + kernels->numStamps); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation", + kernels->mean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", + kernels->rms); } Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionAnalysis.h (revision 25407) @@ -27,4 +27,5 @@ bool pmSubtractionAnalysis( psMetadata *analysis, ///< Metadata container for QA information + psMetadata *header, ///< Metadata container for QA information to put in header pmSubtractionKernels *kernels, ///< Kernels psRegion *region, ///< Region for subtraction Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.c (revision 25407) @@ -765,2 +765,33 @@ return PM_SUBTRACTION_KERNEL_NONE; } + +pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(in, NULL); + + pmSubtractionKernels *out = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return + psMemSetDeallocator(out, (psFreeFunc)subtractionKernelsFree); + + out->type = in->type; + out->description = psMemIncrRefCounter(in->description); + out->num = in->num; + out->u = psMemIncrRefCounter(in->u); + out->v = psMemIncrRefCounter(in->v); + out->widths = psMemIncrRefCounter(in->widths); + out->preCalc = psMemIncrRefCounter(in->preCalc); + out->penalty = in->penalty; + out->penalties = psMemIncrRefCounter(in->penalties); + out->uStop = psMemIncrRefCounter(in->uStop); + out->vStop = psMemIncrRefCounter(in->vStop); + out->size = in->size; + out->inner = in->inner; + out->spatialOrder = in->spatialOrder; + out->bgOrder = in->bgOrder; + out->mode = in->mode; + out->numCols = in->numCols; + out->numRows = in->numRows; + out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL; + out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL; + + return out; +} Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionKernels.h (revision 25407) @@ -203,4 +203,11 @@ ); +/// Copy kernels +/// +/// A deep copy is performed on the solution only; the other components are merely pointers. +pmSubtractionKernels *pmSubtractionKernelsCopy( + const pmSubtractionKernels *in // Kernels to copy + ); + #endif Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c (revision 25407) @@ -10,4 +10,5 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmHDUUtils.h" #include "pmSubtractionParams.h" #include "pmSubtractionKernels.h" @@ -58,16 +59,16 @@ -static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *ro1, // Readout 1 - const pmReadout *ro2, // Readout 2 - const psImage *subMask, // Mask for subtraction, or NULL - psImage *variance, // Variance map - const psRegion *region, // Region of interest, or NULL - float thresh1, // Threshold for stamp finding on readout 1 - float thresh2, // Threshold for stamp finding on readout 2 - float stampSpacing, // Spacing between stamps - int size, // Kernel half-size - int footprint, // Convolution footprint for stamps - pmSubtractionMode mode // Mode for subtraction +static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *variance, // Variance map + const psRegion *region, // Region of interest, or NULL + float thresh1, // Threshold for stamp finding on readout 1 + float thresh2, // Threshold for stamp finding on readout 2 + float stampSpacing, // Spacing between stamps + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction ) { @@ -163,4 +164,31 @@ } +static void subtractionAnalysisUpdate(pmReadout *conv1, pmReadout *conv2, // Convolved images + const psMetadata *analysis, // Analysis metadata + const psMetadata *header // Header metadata + ) +{ + if (conv1) { + conv1->analysis = psMetadataCopy(conv1->analysis, analysis); + } + if (conv2) { + conv2->analysis = psMetadataCopy(conv2->analysis, analysis); + } + + if (conv1 && conv1->parent) { + pmHDU *hdu = pmHDUFromCell(conv1->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, header); + } + } + if (conv2 && conv2->parent) { + pmHDU *hdu = pmHDUFromCell(conv2->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, header); + } + } + + return; +} @@ -253,4 +281,5 @@ psMetadata *outAnalysis = psMetadataAlloc(); // Output analysis values + psMetadata *outHeader = psMetadataAlloc(); // Output header values psTrace("psModules.imcombine", 2, "Convolving...\n"); @@ -259,7 +288,9 @@ psRegion *region = regions->data[i]; // Region of interest - if (!pmSubtractionAnalysis(outAnalysis, kernel, region, ro1->image->numCols, ro1->image->numRows)) { + if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, + ro1->image->numCols, ro1->image->numRows)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data"); psFree(outAnalysis); + psFree(outHeader); psFree(subMask); psFree(kernels); @@ -272,4 +303,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); psFree(outAnalysis); + psFree(outHeader); psFree(subMask); psFree(kernels); @@ -283,11 +315,7 @@ psFree(regions); - if (conv1) { - psMetadataCopy(conv1->analysis, outAnalysis); - } - if (conv2) { - psMetadataCopy(conv2->analysis, outAnalysis); - } + subtractionAnalysisUpdate(conv1, conv2, outAnalysis, outHeader); psFree(outAnalysis); + psFree(outHeader); return true; @@ -369,4 +397,5 @@ pmSubtractionKernels *kernels = NULL; // Kernel basis functions psMetadata *analysis = psMetadataAlloc(); // QA data + psMetadata *header = psMetadataAlloc(); // QA data for header int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions @@ -442,10 +471,29 @@ // We get the stamps here; we will also attempt to get stamps at the first iteration, but it // doesn't matter. - if (!getStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, + if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, stampSpacing, size, footprint, subMode)) { goto MATCH_ERROR; } + + // Define kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, + stamps, footprint, optThreshold, penalty, subMode); + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder, penalty, subMode); + } + + memCheck("kernels"); + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { +#if 0 // Get backgrounds psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background @@ -469,4 +517,6 @@ pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use +#endif + pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { case PM_SUBTRACTION_MODE_1: @@ -483,27 +533,11 @@ } - // Define kernel basis functions - if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold, penalty, subMode); - if (!kernels) { - psErrorClear(); - psWarning("Unable to derive optimum ISIS kernel --- switching to default."); - } - } - if (kernels == NULL) { - // Not an ISIS/GUNK kernel, or the optimum kernel search failed - kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, subMode); - } - - memCheck("kernels"); - int numRejected = -1; // Number of rejected stamps in each iteration for (int k = 0; k < iter && numRejected != 0; k++) { psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, - stampSpacing, size, footprint, subMode)) { + if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, + stampThresh1, stampThresh2, stampSpacing, + size, footprint, subMode)) { goto MATCH_ERROR; } @@ -535,5 +569,5 @@ psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint); + numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); @@ -557,5 +591,5 @@ goto MATCH_ERROR; } - pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint); + pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); psFree(deviations); } @@ -565,5 +599,5 @@ memCheck("solution"); - if (!pmSubtractionAnalysis(analysis, kernels, region, numCols, numRows)) { + if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data"); goto MATCH_ERROR; @@ -601,11 +635,7 @@ memCheck("convolution"); - if (conv1) { - psMetadataCopy(conv1->analysis, analysis); - } - if (conv2) { - psMetadataCopy(conv2->analysis, analysis); - } + subtractionAnalysisUpdate(conv1, conv2, analysis, header); psFree(analysis); + psFree(header); #ifdef TESTING @@ -629,4 +659,5 @@ MATCH_ERROR: psFree(analysis); + psFree(header); psFree(region); psFree(regionString); @@ -842,2 +873,137 @@ return mode; } + + +// Test a subtraction mode by performing a single iteration +static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode + pmSubtractionKernels *kernels, // Kernel description + const char *description, // Description for trace + psImage *subMask, // Subtraction mask + float rej // Rejection threshold + ) +{ + assert(stamps); + assert(kernels); + + psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return false; + } + + psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); + + if (numRejected > 0) { + // Allow re-fit with reduced stamps set + psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return false; + } + psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); + } + + return true; +} + + +pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels, + const psImage *subMask, float rej) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR); + + // Copies of the inputs + pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type); + kernels1->mode = PM_SUBTRACTION_MODE_1; + + if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 1"); + psFree(stamps1); + psFree(kernels1); + psFree(subMask1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask1); + + // Copies of the inputs + pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type); + kernels2->mode = PM_SUBTRACTION_MODE_2; + + if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test subtraction with convolution of image 2"); + psFree(stamps2); + psFree(kernels2); + psFree(subMask2); + psFree(stamps1); + psFree(kernels1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask2); + + + pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps + pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n", + kernels1->mean, kernels1->rms, kernels1->numStamps, + kernels2->mean, kernels2->rms, kernels2->numStamps); + + if (kernels1->mean < kernels2->mean) { + bestStamps = stamps1; + bestKernels = kernels1; + } else { + bestStamps = stamps2; + bestKernels = kernels2; + } + + psFree(*stamps); + psFree(*kernels); + *stamps = psMemIncrRefCounter(bestStamps); + *kernels = psMemIncrRefCounter(bestKernels); + + psFree(stamps1); + psFree(stamps2); + psFree(kernels1); + psFree(kernels2); + + return bestKernels->mode; +} + Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.h (revision 25407) @@ -83,4 +83,13 @@ ); +/// Determine best subtraction mode to use +/// +/// Subtractions are attempted each way, and the mode with the lower residual is taken to be the best +pmSubtractionMode pmSubtractionBestMode( + pmSubtractionStampList **stamps, ///< Stamps to use for solution + pmSubtractionKernels **kernels, ///< Kernels to use for solution + const psImage *subMask, ///< Subtraction mask + float rej ///< Rejection threshold for stamps + ); #endif Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.c (revision 25407) @@ -225,4 +225,70 @@ } +pmSubtractionStampList *pmSubtractionStampListCopy(const pmSubtractionStampList *in) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(in, NULL); + + pmSubtractionStampList *out = psAlloc(sizeof(pmSubtractionStampList)); // Copied stamp list to return + psMemSetDeallocator(out, (psFreeFunc)subtractionStampListFree); + + int num = out->num = in->num; // Number of stamps + out->stamps = psArrayAlloc(num); + out->regions = psArrayAlloc(num); + out->x = NULL; + out->y = NULL; + out->flux = NULL; + out->footprint = in->footprint; + + for (int i = 0; i < num; i++) { + psRegion *inRegion = in->regions->data[i]; // Input region + out->regions->data[i] = psRegionAlloc(inRegion->x0, inRegion->x1, inRegion->y0, inRegion->y1); + + pmSubtractionStamp *inStamp = in->stamps->data[i]; // Input stamp + pmSubtractionStamp *outStamp = psAlloc(sizeof(pmSubtractionStamp)); + psMemSetDeallocator(outStamp, (psFreeFunc)subtractionStampFree); + outStamp->x = inStamp->x; + outStamp->y = inStamp->y; + outStamp->flux = inStamp->flux; + outStamp->xNorm = inStamp->xNorm; + outStamp->yNorm = inStamp->yNorm; + outStamp->status = inStamp->status; + + outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL; + outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL; + outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL; + + if (inStamp->convolutions1) { + int size = inStamp->convolutions1->n; // Size of array + outStamp->convolutions1 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + psKernel *conv = inStamp->convolutions1->data[j]; // Convolution + outStamp->convolutions1->data[j] = conv ? psKernelCopy(conv) : NULL; + } + } else { + outStamp->convolutions1 = NULL; + } + if (inStamp->convolutions2) { + int size = inStamp->convolutions2->n; // Size of array + outStamp->convolutions2 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + psKernel *conv = inStamp->convolutions2->data[j]; // Convolution + outStamp->convolutions2->data[j] = conv ? psKernelCopy(conv) : NULL; + } + } else { + outStamp->convolutions2 = NULL; + } + + outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL; + outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL; + outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL; + outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL; + outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL; + + out->stamps->data[i] = outStamp; + } + + return out; +} + pmSubtractionStamp *pmSubtractionStampAlloc(void) { Index: /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionStamps.h (revision 25407) @@ -52,4 +52,12 @@ } \ } + +/// Copy a list of stamps +/// +/// A deep copy is performed of the stamp list and the component stamps +pmSubtractionStampList *pmSubtractionStampListCopy( + const pmSubtractionStampList *in // Stamp list to copy + ); + /// A stamp for image subtraction Index: /branches/eam_branches/20090715/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/Makefile.am (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/objects/Makefile.am (revision 25407) @@ -53,5 +53,6 @@ pmGrowthCurveGenerate.c \ pmGrowthCurve.c \ - pmSourceMatch.c + pmSourceMatch.c \ + pmDetEff.c EXTRA_DIST = \ @@ -90,5 +91,6 @@ pmTrend2D.h \ pmGrowthCurve.h \ - pmSourceMatch.h + pmSourceMatch.h \ + pmDetEff.h CLEANFILES = *~ Index: /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.c (revision 25407) +++ /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.c (revision 25407) @@ -0,0 +1,168 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#include "pmDetEff.h" + + +static void detEffFree(pmDetEff *de) +{ + psFree(de->magOffsets); + psFree(de->counts); + psFree(de->magDiffMean); + psFree(de->magDiffStdev); + psFree(de->magErrMean); +} + + +pmDetEff *pmDetEffAlloc(float magRef, int numSources, int numBins) +{ + pmDetEff *de = psAlloc(sizeof(pmDetEff)); // Detection efficiency, to return + psMemSetDeallocator(de, (psFreeFunc)detEffFree); + + de->magRef = magRef; + de->numSources = numSources; + de->numBins = numBins; + + de->magOffsets = NULL; + de->counts = NULL; + de->magDiffMean = NULL; + de->magDiffStdev = NULL; + de->magErrMean = NULL; + + return de; +} + + +bool pmDetEffWrite(psFits *fits, pmDetEff *de, const psMetadata *header, const char *extname) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_FITS_WRITABLE(fits, false); + PM_ASSERT_DETEFF_RESULTS(de, false); + + psArray *table = psArrayAlloc(de->numBins); // Table to write + for (int i = 0; i < de->numBins; i++) { + psMetadata *row = psMetadataAlloc(); // Table row + psMetadataAddF32(row, PS_LIST_TAIL, "OFFSET", 0, "Magnitude offset from reference", + de->magOffsets->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "COUNTS", 0, "Number of sources recovered", + de->counts->data.S32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.MEAN", 0, "Mean magnitude difference", + de->magDiffMean->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.STDEV", 0, "Stdev magnitude difference", + de->magDiffStdev->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ERR.MEAN", 0, "Mean magnitude error", + de->magErrMean->data.F32[i]); + table->data[i] = row; + } + + psMetadata *deHeader = psMetadataCopy(NULL, header); // Header for detection efficiency + psMetadataAddF32(deHeader, PS_LIST_TAIL, "DETEFF.MAGREF", PS_META_REPLACE, "Magnitude reference", + de->magRef); + psMetadataAddS32(deHeader, PS_LIST_TAIL, "DETEFF.NUM", PS_META_REPLACE, "Number of fake sources injected", + de->numSources); + + if (!psFitsWriteTable(fits, deHeader, table, extname)) { + psError(PS_ERR_IO, false, "Unable to write detection efficiency table."); + psFree(table); + psFree(deHeader); + return false; + } + + psFree(table); + psFree(deHeader); + + return true; +} + +bool pmReadoutWriteDetEff(psFits *fits, const pmReadout *readout, + const psMetadata *header, const char *extname) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + + bool mdok; // Status of MD lookup + pmDetEff *de = psMetadataLookupPtr(&mdok, readout->analysis, PM_DETEFF_ANALYSIS); // Detection efficiency + if (!mdok || !de) { + // Wrote everything there was to write + return true; + } + return pmDetEffWrite(fits, de, header, extname); +} + + +pmDetEff *pmDetEffRead(psFits *fits, const char *extname) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_STRING_NON_EMPTY(extname, false); + + if (!psFitsMoveExtNameClean(fits, extname)) { + // Nothing to read + return NULL; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // Header for table + if (!header) { + psError(PS_ERR_IO, false, "Unable to read FITS header"); + return NULL; + } + + int numBins = psFitsTableSize(fits);// Size of table + bool mdok; // Status of MD lookup + int numSources = psMetadataLookupS32(&mdok, header, "DETEFF.NUM"); // Number of fake sources injected + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find number of sources"); + psFree(header); + return NULL; + } + float magRef = psMetadataLookupF32(&mdok, header, "DETEFF.MAGREF"); // Magnitude reference + if (!mdok) { + psError(PS_ERR_UNKNOWN, false, "Unable to find magnitude reference"); + psFree(header); + return NULL; + } + psFree(header); + + pmDetEff *de = pmDetEffAlloc(magRef, numSources, numBins); // Detection efficiency + de->magOffsets = psVectorAlloc(numBins, PS_TYPE_F32); + de->counts = psVectorAlloc(numBins, PS_TYPE_S32); + de->magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); + de->magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); + de->magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); + + psArray *table = psFitsReadTable(fits); // FITS table + if (!table) { + psError(PS_ERR_IO, false, "Unable to read detection efficiency table."); + psFree(de); + return false; + } + + for (int i = 0; i < numBins; i++) { + psMetadata *row = table->data[i]; // Table row + de->magOffsets->data.F32[i] = psMetadataLookupF32(NULL, row, "OFFSET"); + de->counts->data.S32[i] = psMetadataLookupS32(NULL, row, "COUNTS"); + de->magDiffMean->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.MEAN"); + de->magDiffStdev->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.STDEV"); + de->magErrMean->data.F32[i] = psMetadataLookupF32(NULL, row, "ERR.MEAN"); + } + + psFree(table); + return de; +} + +bool pmReadoutReadDetEff(psFits *fits, const pmReadout *readout, const char *extname) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + + pmDetEff *de = pmDetEffRead(fits, extname); + if (!de) { + if (psErrorCodeLast() != PS_ERR_NONE) { + return false; + } + return true; + } + + return psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PM_DETEFF_ANALYSIS, + PS_META_REPLACE | PS_DATA_UNKNOWN, "Detection efficiency", de); +} Index: /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.h (revision 25407) +++ /branches/eam_branches/20090715/psModules/src/objects/pmDetEff.h (revision 25407) @@ -0,0 +1,81 @@ +#ifndef PM_DET_EFF_H +#define PM_DET_EFF_H + +#include +#include + +#include "pmFPA.h" + +#define PM_DETEFF_ANALYSIS "DETEFF" // Location of detection efficiency on pmReadout.analysis + +// Detection efficiency characterisation +typedef struct { + float magRef; // Reference magnitude + int numSources; // Number of sources + int numBins; // Number of bins + psVector *magOffsets; // Magnitude offsets for each bin + psVector *counts; // Counts of sources recovered for each bin + psVector *magDiffMean; // Mean magnitude difference for each bin + psVector *magDiffStdev; // Stdev of magnitude difference for each bin + psVector *magErrMean; // Mean magnitude error for each bin +} pmDetEff; + + +/// Allocator +pmDetEff *pmDetEffAlloc(float magRef, // Reference magnitude + int numSources, // Number of sources + int numBins // Number of bins + ); + +/// Write detection efficiency to FITS file +bool pmDetEffWrite(psFits *fits, // FITS file to which to write + pmDetEff *deteff, // Detection efficiency to write + const psMetadata *header, // Header to write + const char *extname // Extension name + ); + +/// Write detection efficiency from a readout to a FITS file +bool pmReadoutWriteDetEff(psFits *fits,// FITS file to which to write + const pmReadout *readout, // Readout with detection efficiency + const psMetadata *header, // Header to write + const char *extname // Extension name + ); + +/// Read detection efficiency +pmDetEff *pmDetEffRead(psFits *fits, // FITS file from which to read + const char *extname // Extension name + ); + +/// Read detection efficiency into a readout +bool pmReadoutReadDetEff(psFits *fits,// FITS file to which to write + const pmReadout *readout, // Readout with detection efficiency + const char *extname // Extension name + ); + +#define PM_ASSERT_DETEFF_NON_NULL(DE, RETURN) { \ + if (!(DE)) { \ + psError(PS_ERR_UNEXPECTED_NULL, true, "Detection efficiency %s is NULL", #DE); \ + return RETURN; \ + } \ +} + +#define PM_ASSERT_DETEFF_RESULTS(DE, RETURN) { \ + PM_ASSERT_DETEFF_NON_NULL(DE, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magOffsets, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magOffsets, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magOffsets, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->counts, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->counts, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->counts, PS_TYPE_S32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffMean, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magDiffMean, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magDiffMean, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffStdev, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magDiffStdev, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magDiffStdev, PS_TYPE_F32, RETURN); \ + PS_ASSERT_VECTOR_NON_NULL((DE)->magErrMean, RETURN); \ + PS_ASSERT_VECTOR_SIZE((DE)->magErrMean, (long)(DE)->numBins, RETURN); \ + PS_ASSERT_VECTOR_TYPE((DE)->magErrMean, PS_TYPE_F32, RETURN); \ +} + +#endif Index: /branches/eam_branches/20090715/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmSourceIO.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/objects/pmSourceIO.c (revision 25407) @@ -42,8 +42,82 @@ #include "pmSource.h" #include "pmModelClass.h" +#include "pmDetEff.h" #include "pmSourceIO.h" #define BLANK_HEADERS "BLANK.HEADERS" // Name of metadata in camera configuration containing header names // for putting values into a blank PHU + +// lookup the EXTNAME values used for table data and image header segments +static bool sourceExtensions(psString *headname, // Extension name for header + psString *dataname, // Extension name for data + psString *deteffname, // Extension name for detection efficiency + psString *xsrcname, // Extension name for extended sources + psString *xfitname, // Extension name for extended fits + const pmFPAfile *file, // File of interest + const pmFPAview *view // View to level of interest + ) +{ + bool status; // Status of MD lookup + + // Menu of EXTNAME rules + psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); + if (!menu) { + psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); + return false; + } + + // EXTNAME for image header + if (headname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); + return false; + } + *headname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for table data + if (dataname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); + return false; + } + *dataname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for detection efficiency + if (deteffname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.DETEFF"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DETEFF in EXTNAME.RULES in camera.config"); + return false; + } + *deteffname = pmFPAfileNameFromRule(rule, file, view); + } + + // EXTNAME for extended source data table + if (xsrcname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config"); + return false; + } + *xsrcname = pmFPAfileNameFromRule (rule, file, view); + } + + if (xfitname) { + // EXTNAME for extended source data table + const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config"); + return false; + } + *xfitname = pmFPAfileNameFromRule (rule, file, view); + } + + return true; +} + // translations between psphot object types and dophot object types @@ -271,8 +345,4 @@ char *exttype = NULL; - char *dataname = NULL; - char *xsrcname = NULL; - char *xfitname = NULL; - char *headname = NULL; // if sources is NULL, write out an empty table @@ -354,49 +424,12 @@ // define the EXTNAME values for the different data segments: - { - // lookup the EXTNAME values used for table data and image header segments - char *rule = NULL; - - // Menu of EXTNAME rules - psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); - if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; - } - - // EXTNAME for image header - rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); - return false; - } - headname = pmFPAfileNameFromRule (rule, file, view); - - // EXTNAME for table data - rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); - return false; - } - dataname = pmFPAfileNameFromRule (rule, file, view); - - if (XSRC_OUTPUT) { - // EXTNAME for extended source data table - rule = psMetadataLookupStr(&status, menu, "CMF.XSRC"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config"); - return false; - } - xsrcname = pmFPAfileNameFromRule (rule, file, view); - } - if (XFIT_OUTPUT) { - // EXTNAME for extended source data table - rule = psMetadataLookupStr(&status, menu, "CMF.XFIT"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config"); - return false; - } - xfitname = pmFPAfileNameFromRule (rule, file, view); - } + psString headname = NULL; + psString dataname = NULL; + psString deteffname = NULL; + psString xsrcname = NULL; + psString xfitname = NULL; + if (!sourceExtensions(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL, + XFIT_OUTPUT ? &xfitname : NULL, file, view)) { + return false; } @@ -480,49 +513,54 @@ // XXX these are case-sensitive since the EXTYPE is case-sensitive - status = false; + status = true; if (!strcmp (exttype, "SMPDATA")) { - status = pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_DEV_0")) { - status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname); + status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); + 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); - } + status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); + } + + if (deteffname) { + status &= pmReadoutWriteDetEff(file->fits, readout, outhead, deteffname); + } + if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); + status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); + 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); + status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); } } if (xfitname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname); } if (!strcmp (exttype, "PS1_CAL_0")) { - status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname); + status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname); } if (!strcmp (exttype, "PS1_V1")) { - status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); + 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); + status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); } } @@ -572,6 +610,6 @@ // not needed if only one chip if (file->fpa->chips->n == 1) { - pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); - return true; + pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config); + return true; } @@ -885,26 +923,11 @@ hdu = pmFPAviewThisHDU (view, file->fpa); - // lookup the EXTNAME values used for table data and image header segments - char *rule = NULL; - // Menu of EXTNAME rules - psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); - if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; - } - // EXTNAME for image header - rule = psMetadataLookupStr(&status, menu, "CMF.HEAD"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config"); - return false; - } - char *headname = pmFPAfileNameFromRule (rule, file, view); - // EXTNAME for table data - rule = psMetadataLookupStr(&status, menu, "CMF.DATA"); - if (!rule) { - psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config"); - return false; - } - char *dataname = pmFPAfileNameFromRule (rule, file, view); + // define the EXTNAME values for the different data segments: + psString headname = NULL; + psString dataname = NULL; + psString deteffname = NULL; + if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) { + return false; + } // advance to the IMAGE HEADER extension @@ -958,4 +981,9 @@ sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); } + + if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) { + psError(PS_ERR_IO, false, "Unable to read detection efficiency"); + return false; + } } @@ -1070,3 +1098,3 @@ } - + Index: /branches/eam_branches/20090715/psModules/src/objects/pmSourceMatch.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmSourceMatch.c (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/objects/pmSourceMatch.c (revision 25407) @@ -221,5 +221,5 @@ } else { // Match with the master list - psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources + psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree long numMatch = 0; // Number of matches Index: /branches/eam_branches/20090715/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/psmodules.h (revision 25406) +++ /branches/eam_branches/20090715/psModules/src/psmodules.h (revision 25407) @@ -133,4 +133,5 @@ #include #include +#include // The following headers are from random locations, here because they cross bounds