Index: branches/eam_branches/20090522/psModules/src/camera/pmFPABin.c =================================================================== --- branches/eam_branches/20090522/psModules/src/camera/pmFPABin.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/camera/pmFPABin.c (revision 24557) @@ -65,5 +65,7 @@ } - float imageValue, maskValue;// Values to set + // Values to set + float imageValue; + psImageMaskType maskValue; if (numPix > 0) { imageValue = sum / numPix; Index: branches/eam_branches/20090522/psModules/src/camera/pmFPACopy.c =================================================================== --- branches/eam_branches/20090522/psModules/src/camera/pmFPACopy.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/camera/pmFPACopy.c (revision 24557) @@ -70,4 +70,5 @@ psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); *target = psImageAlloc(binning->nXruff, binning->nYruff, source->type.type); + psImageInit (*target, 0.0); return; } @@ -221,4 +222,11 @@ psTrace("psModules.camera", 3, "xFlip: %d; yFlip: %d\n", xFlip, yFlip); + // Blow away extant readouts + for (int i = 0; i < target->readouts->n; i++) { + psFree(target->readouts->data[i]); + target->readouts->data[i] = NULL; + } + target->readouts->n = 0; + // Perform deep copy of the images. I would prefer *not* to do a deep copy, in the interests of speed (we // still need to do another deep copy into the HDU for when we write out), but this is the only way I can @@ -242,4 +250,21 @@ readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip, pixels); + // Copy covariance matrix: doesn't care about flips, etc. + if (sourceReadout->covariance) { + if (targetReadout->covariance) { + psFree(targetReadout->covariance); + } + targetReadout->covariance = psKernelCopy(sourceReadout->covariance); +#if 0 + if (binning) { + // XXX This isn't strictly correct, but we don't have a function that bins covariance matrices + // with unequal binning factors. + psKernel *covar = psImageCovarianceBin(PS_MAX(binning->nXbin, binning->nYbin), + targetReadout->covariance); + psFree(targetReadout->covariance); + targetReadout->covariance = covar; + } +#endif + } // Copy bias Index: branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.c (revision 24557) @@ -199,7 +199,8 @@ } -bool pmReadoutSetVariance(pmReadout *readout, bool poisson) +bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); + // check that the noiseMap (if it exists) matches the readout variance size) pmCell *cell = readout->parent; // The parent cell @@ -228,4 +229,5 @@ // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1. + // XXX this calculation is wrong: limit is 1 e-, but this is in DN readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs"); readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max", @@ -239,6 +241,12 @@ } - readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", - psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + // apply a supplied readnoise map (NOTE: in DN, not electrons): + if (noiseMap) { + psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap); + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar); + psFree (rdVar); + } else { + readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32)); + } return true; @@ -247,5 +255,5 @@ // this function creates the variance pixels, or uses the existing variance pixels. it will set // the noise pixel values only if the variance image is not supplied -bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson) +bool pmReadoutGenerateVariance(pmReadout *readout, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -291,8 +299,8 @@ readout->variance = variance; - return pmReadoutSetVariance(readout, poisson); -} - -bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson) + return pmReadoutSetVariance(readout, noiseMap, poisson); +} + +bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -301,10 +309,10 @@ success &= pmReadoutGenerateMask(readout, satMask, badMask); - success &= pmReadoutGenerateVariance(readout, poisson); + success &= pmReadoutGenerateVariance(readout, noiseMap, poisson); return success; } -bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson) +bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson) { PS_ASSERT_PTR_NON_NULL(cell, false); @@ -314,5 +322,5 @@ for (int i = 0; i < readouts->n; i++) { pmReadout *readout = readouts->data[i]; // The readout - success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask); + success &= pmReadoutGenerateMaskVariance(readout, satMask, badMask, noiseMap, poisson); } Index: branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.h (revision 24529) +++ branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.h (revision 24557) @@ -54,4 +54,5 @@ /// can't be generated. bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -72,4 +73,5 @@ /// with HDU entry). This is intended for most operations. bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -81,4 +83,5 @@ psImageMaskType sat, ///< Mask value to give saturated pixels psImageMaskType bad, ///< Mask value to give bad (low) pixels + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); @@ -90,4 +93,5 @@ psImageMaskType sat, ///< Mask value to give saturated pixels psImageMaskType bad, ///< Mask value to give bad (low) pixels + const psImage *noiseMap, ///< 2D image of the read noise in DN bool poisson ///< Include poisson variance (in addition to read noise)? ); Index: branches/eam_branches/20090522/psModules/src/camera/pmHDU.c =================================================================== --- branches/eam_branches/20090522/psModules/src/camera/pmHDU.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/camera/pmHDU.c (revision 24557) @@ -93,12 +93,31 @@ } + psTrace("psModules.camera", 5, "Reading the header...\n"); + + // The header may already exist (e.g., from doing concept writing at the PHU level) so we need to be + // careful. We read into a separate container and copy that over the top of anything that's already read. + psMetadata *header = psFitsReadHeader(NULL, fits); + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); + return false; + } + if (!hdu->header) { - psTrace("psModules.camera", 5, "Reading the header...\n"); - hdu->header = psFitsReadHeader(hdu->header, fits); - if (! hdu->header) { - psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); - return false; - } - } + hdu->header = header; + return true; + } + + psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + const char *name = item->name; // Name of item + if (psMetadataLookup(hdu->header, name)) { + // It exists; clobber + psMetadataRemoveKey(hdu->header, name); + } + psMetadataAddItem(hdu->header, item, PS_LIST_TAIL, 0); + } + psFree(iter); + psFree(header); return true; Index: branches/eam_branches/20090522/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- branches/eam_branches/20090522/psModules/src/concepts/pmConceptsStandard.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/concepts/pmConceptsStandard.c (revision 24557) @@ -36,9 +36,9 @@ // Format type for time typedef enum { - TIME_FORMAT_YYYYMMDD, // Date stored in YYYY-MM-DD order (ISO-standard) - TIME_FORMAT_DDMMYYYY, // Date stored in DD-MM-YYYY order - TIME_FORMAT_MMDDYYYY, // Date stored in MM-DD-YYYY order - TIME_FORMAT_JD, // Date stored as JD - TIME_FORMAT_MJD, // Date stored as MJD + TIME_FORMAT_YYYYMMDD, // Date stored in YYYY-MM-DD order (ISO-standard) + TIME_FORMAT_DDMMYYYY, // Date stored in DD-MM-YYYY order + TIME_FORMAT_MMDDYYYY, // Date stored in MM-DD-YYYY order + TIME_FORMAT_JD, // Date stored as JD + TIME_FORMAT_MJD, // Date stored as MJD } conceptTimeFormatType; @@ -391,4 +391,5 @@ } + // FPA.RA and FPA.DEC psMetadataItem *p_pmConceptFormat_FPA_Coords(const psMetadataItem *concept, @@ -410,9 +411,8 @@ // How to interpret the coordinates bool mdok = true; // Status of MD lookup - psMetadata *formats = psMetadataLookupMetadata(&mdok, - cameraFormat, - "FORMATS"); + psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS"); + bool sexagesimal = false; // Write sexagesimal format? if (mdok && formats) { - psString format = psMetadataLookupStr(&mdok,formats, concept->name); + psString format = psMetadataLookupStr(&mdok, formats, concept->name); if (mdok && strlen(format) > 0) { if (strcasecmp(format, "HOURS") == 0) { @@ -428,26 +428,36 @@ coords /= defaultCoordScaling(concept); } + + psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA + psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec + if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) { + sexagesimal = true; + } } else { coords /= defaultCoordScaling(concept); } - // We choose to write sexagesimal format - int big, medium; // Degrees and minutes - float small; // Seconds - bool negative = (coords < 0); // Are we working below zero? - coords = fabs(coords); - big = (int)abs(coords); - coords -= big; - medium = 60.0 * coords; - coords -= medium / 60.0; - small = 3600.0 * coords; - small = (float)((int)(small * 1000.0)) / 1000.0; - if (negative) { - big *= -1; - } - psString coordString = NULL; // String with the coordinates in sexagesimal format - psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small); - psMetadataItem *coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString); - psFree(coordString); + psMetadataItem *coordItem = NULL; // Item with coordinates, to return + if (sexagesimal) { + int big, medium; // Degrees and minutes + float small; // Seconds + bool negative = (coords < 0); // Are we working below zero? + coords = fabs(coords); + big = (int)abs(coords); + coords -= big; + medium = 60.0 * coords; + coords -= medium / 60.0; + small = 3600.0 * coords; + small = (float)((int)(small * 1000.0)) / 1000.0; + if (negative) { + big *= -1; + } + psString coordString = NULL; // String with the coordinates in sexagesimal format + psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small); + coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString); + psFree(coordString); + } else { + coordItem = psMetadataItemAllocF64(concept->name, concept->comment, coords); + } return coordItem; Index: branches/eam_branches/20090522/psModules/src/config/pmConfig.c =================================================================== --- branches/eam_branches/20090522/psModules/src/config/pmConfig.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/config/pmConfig.c (revision 24557) @@ -243,4 +243,6 @@ unsigned int numBadLines = 0; struct stat filestat; + + pmErrorRegister(); psTrace("psModules.config", 3, "Loading %s configuration from file %s\n", @@ -1496,4 +1498,5 @@ } } + globfree(&globList); } psFree (words); Index: branches/eam_branches/20090522/psModules/src/config/pmErrorCodes.dat =================================================================== --- branches/eam_branches/20090522/psModules/src/config/pmErrorCodes.dat (revision 24529) +++ branches/eam_branches/20090522/psModules/src/config/pmErrorCodes.dat (revision 24557) @@ -12,4 +12,5 @@ SKY Problem in sky STAMPS Unable to select stamps for PSF-matching +SMALL_AREA Small area for PSF-matching # these errors correspond to standard exit conditions ARGUMENTS Incorrect arguments Index: branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.c =================================================================== --- branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.c (revision 24557) @@ -105,4 +105,5 @@ DETREND_STRING_CASE(FRINGE); DETREND_STRING_CASE(ASTROM); + DETREND_STRING_CASE(NOISEMAP); default: return NULL; Index: branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.h =================================================================== --- branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.h (revision 24529) +++ branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.h (revision 24557) @@ -35,4 +35,5 @@ PM_DETREND_TYPE_BACKGROUND, PM_DETREND_TYPE_ASTROM, + PM_DETREND_TYPE_NOISEMAP, } pmDetrendType; Index: branches/eam_branches/20090522/psModules/src/detrend/pmDetrendThreads.c =================================================================== --- branches/eam_branches/20090522/psModules/src/detrend/pmDetrendThreads.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/detrend/pmDetrendThreads.c (revision 24557) @@ -60,4 +60,6 @@ } + // NOISEMAP : for now, not applied in the threaded loop + return true; } Index: branches/eam_branches/20090522/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmSubtraction.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmSubtraction.c (revision 24557) @@ -196,49 +196,12 @@ { psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image - int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy for (int y = -footprint; y <= footprint; y++) { - // Convolution with a delta function is just the value specified by the offset - memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes); + for (int x = -footprint; x <= footprint; x++) { + convolved->kernel[y][x] = image->kernel[y - v][x - u]; + } } return convolved; } -// Take a subset of an image -static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved - psRegion region, // Region of interest (with border) - bool threaded // Are we running threaded? - ) -{ - if (!image) { - return NULL; - } - // XXX if (threaded) { - // XXX psMutexLock(image); - // XXX } - psImage *subset = psImageSubset(image, region); // Subset image, to return - // XXX if (threaded) { - // XXX psMutexUnlock(image); - // XXX } - return subset; -} - -// Free the subset of an image -static inline void convolveSubsetFree(psImage *parent, // Parent image - psImage *child, // Child (subset) image - bool threaded // Are we running threaded? - ) -{ - if (!child || !parent) { - return; - } - // XXX if (threaded) { - // XXX psMutexLock(parent); - // XXX } - psFree(child); - // XXX if (threaded) { - // XXX psMutexUnlock(parent); - // XXX } - return; -} // Convolve an image using FFT @@ -256,13 +219,11 @@ region.y0 - size, region.y1 + size); // Add a border - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - - psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve - psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask + psImage *subImage = image ? psImageSubset(image, border) : NULL; // Subimage to convolve + psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution - convolveSubsetFree(image, subImage, threaded); - convolveSubsetFree(mask, subMask, threaded); + psFree(subImage); + psFree(subMask); // Now, we have to stick it in where it belongs @@ -300,9 +261,7 @@ region.y0 - size, region.y1 + size); // Add a border - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - - psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map - psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image - psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask + psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map + psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image + psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once @@ -310,7 +269,7 @@ psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys - convolveSubsetFree(variance, subVariance, threaded); - convolveSubsetFree(sys, subSys, threaded); - convolveSubsetFree(mask, subMask, threaded); + psFree(subVariance); + psFree(subSys); + psFree(subMask); // Now, we have to stick it in where it belongs @@ -420,15 +379,13 @@ if (box > 0) { int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds - bool threaded = pmSubtractionThreaded(); // Are we running threaded? - psRegion region = psRegionSet(colMin - box, colMax + box, rowMin - box, rowMax + box); // Region to convolve - psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve + psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad, -box, box, -box, box); // Convolved subtraction mask - convolveSubsetFree(subMask, image, threaded); + psFree(image); psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns"); @@ -653,7 +610,7 @@ float value = 0.0; // Value of convolved pixel int uMin = x - size, uMax = x + size; // Range for u - psF32 *xKernelData = xKernel->data.F32; // Kernel values + psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values psF32 *imageData = &image->kernel[y][uMin]; // Image values - for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) { + for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { value += *xKernelData * *imageData; } @@ -668,7 +625,7 @@ float value = 0.0; // Value of convolved pixel int vMin = y - size, vMax = y + size; // Range for v - psF32 *yKernelData = yKernel->data.F32; // Kernel values + psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! - for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) { + for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { value += *yKernelData * *imageData; } Index: branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionEquation.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionEquation.c (revision 24557) @@ -15,5 +15,7 @@ #include "pmSubtractionVisual.h" -//#define TESTING +//#define TESTING // TESTING output for debugging; may not work with threads! + +#define USE_VARIANCE // Include variance in equation? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -31,5 +33,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = image1->kernel[y][x] * image2->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sum += value; } } @@ -195,5 +201,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumC += value; } } @@ -218,5 +228,8 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; + double invNoise2 = 1.0; +#ifdef USE_VARIANCE + invNoise2 /= variance->kernel[y][x]; +#endif double value = input->kernel[y][x] * invNoise2; sumI += value; @@ -277,5 +290,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = target->kernel[y][x] * conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumTC += value; } } @@ -297,5 +314,8 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = target->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif sumIT += value * input->kernel[y][x]; sumT += value; @@ -366,5 +386,9 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; + double value = conv->kernel[y][x]; +#ifdef USE_VARIANCE + value /= variance->kernel[y][x]; +#endif + sumC += value; } } Index: branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionKernels.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionKernels.c (revision 24557) @@ -145,10 +145,8 @@ // Calculate moments - double sum = 0.0; // Sum of kernel component, for normalisation double moment = 0.0; // Moment, for penalty for (int v = -size, y = 0; v <= size; v++, y++) { for (int u = -size, x = 0; u <= size; u++, x++) { double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel - sum += value; moment += value * (PS_SQR(u) + PS_SQR(v)); } @@ -163,8 +161,8 @@ } } - sum = 1.0 / sum; + sum = 1.0 / sqrt(sum); psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); - moment *= sum; + moment *= PS_SQR(sum); } Index: branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMask.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMask.c (revision 24557) @@ -6,4 +6,5 @@ #include +#include "pmErrorCodes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -78,5 +79,5 @@ } if (numBad > badFrac * numCols * numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, + psError(PM_ERR_SMALL_AREA, true, "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); Index: branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMatch.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMatch.c (revision 24557) @@ -246,5 +246,5 @@ badFrac, mode); // Subtraction mask if (!subMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psError(psErrorCodeLast(), false, "Unable to generate subtraction mask."); psFree(kernels); psFree(regions); @@ -337,5 +337,5 @@ PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false); } - PS_ASSERT_INT_POSITIVE(iter, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); @@ -379,5 +379,5 @@ badFrac, subMode); if (!subMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psError(psErrorCodeLast(), false, "Unable to generate subtraction mask."); goto MATCH_ERROR; } Index: branches/eam_branches/20090522/psModules/src/objects/Makefile.am =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/Makefile.am (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/Makefile.am (revision 24557) @@ -38,4 +38,5 @@ pmSourceIO_PS1_CAL_0.c \ pmSourceIO_CMF_PS1_V1.c \ + pmSourceIO_CMF_PS1_V2.c \ pmSourceIO_MatchedRefs.c \ pmSourcePlots.c \ Index: branches/eam_branches/20090522/psModules/src/objects/pmFootprintArraysMerge.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmFootprintArraysMerge.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmFootprintArraysMerge.c (revision 24557) @@ -28,18 +28,26 @@ */ psArray *pmFootprintArraysMerge(const psArray *footprints1, // one set of footprints - const psArray *footprints2, // the other set - const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2 - assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); - assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); + const psArray *footprints2, // the other set + const int includePeaks // which peaks to set? 0x1 => footprints1, 0x2 => 2 + ) +{ + if (!footprints1 && !footprints2) { + // No footprints in merged array + return psArrayAllocEmpty(0); + } - if (footprints1->n == 0 || footprints2->n == 0) { // nothing to do but put copies on merged - const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1; + assert(!footprints1 || footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); + assert(!footprints2 || footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); - psArray *merged = psArrayAllocEmpty(old->n); - for (int i = 0; i < old->n; i++) { - psArrayAdd(merged, 1, old->data[i]); - } - - return merged; + if (!footprints1 || footprints1->n == 0 || !footprints2 || footprints2->n == 0) { + // nothing to do but put copies on merged + const psArray *old = (!footprints1 || footprints1->n == 0) ? footprints2 : footprints1; + + psArray *merged = psArrayAllocEmpty(old->n); + for (int i = 0; i < old->n; i++) { + psArrayAdd(merged, 1, old->data[i]); + } + + return merged; } /* @@ -48,14 +56,14 @@ */ { - pmFootprint *fp1 = footprints1->data[0]; - pmFootprint *fp2 = footprints2->data[0]; - if (fp1->region.x0 != fp2->region.x0 || - fp1->region.x1 != fp2->region.x1 || - fp1->region.y0 != fp2->region.y0 || - fp1->region.y1 != fp2->region.y1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "The two pmFootprint arrays correspnond to different-sized regions"); - return NULL; - } + pmFootprint *fp1 = footprints1->data[0]; + pmFootprint *fp2 = footprints2->data[0]; + if (fp1->region.x0 != fp2->region.x0 || + fp1->region.x1 != fp2->region.x1 || + fp1->region.y0 != fp2->region.y0 || + fp1->region.y1 != fp2->region.y1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "The two pmFootprint arrays correspnond to different-sized regions"); + return NULL; + } } /* @@ -72,17 +80,17 @@ * Now assign the peaks appropriately. We could do this more efficiently * using idImage (which we just freed), but this is easy and probably fast enough - */ + */ if (includePeaks & 0x1) { - psArray *peaks = pmFootprintArrayToPeaks(footprints1); - pmFootprintsAssignPeaks(merged, peaks); - psFree(peaks); + psArray *peaks = pmFootprintArrayToPeaks(footprints1); + pmFootprintsAssignPeaks(merged, peaks); + psFree(peaks); } if (includePeaks & 0x2) { - psArray *peaks = pmFootprintArrayToPeaks(footprints2); - pmFootprintsAssignPeaks(merged, peaks); - psFree(peaks); + psArray *peaks = pmFootprintArrayToPeaks(footprints2); + pmFootprintsAssignPeaks(merged, peaks); + psFree(peaks); } - + return merged; } Index: branches/eam_branches/20090522/psModules/src/objects/pmMoments.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmMoments.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmMoments.c (revision 24557) @@ -51,4 +51,5 @@ tmp->Sky = 0.0; tmp->nPixels = 0; + tmp->SN = 0; psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); Index: branches/eam_branches/20090522/psModules/src/objects/pmSource.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSource.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSource.c (revision 24557) @@ -480,14 +480,14 @@ if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats"); + psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats"); return (emptyClump); - } + } psfClump.X = stats->clippedMean; psfClump.dX = stats->clippedStdev; if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats"); + psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats"); return (emptyClump); - } + } psfClump.Y = stats->clippedMean; psfClump.dY = stats->clippedStdev; @@ -548,12 +548,12 @@ pmSource *source = (pmSource *) sources->data[i]; - // psf clumps are found for image subregions: - // skip sources not in this region + // psf clumps are found for image subregions: + // skip sources not in this region if (source->peak->x < region->x0) continue; if (source->peak->x >= region->x1) continue; if (source->peak->y < region->y0) continue; - if (source->peak->y >= region->y1) continue; - - // should be set by pmSourceAlloc + if (source->peak->y >= region->y1) continue; + + // should be set by pmSourceAlloc psAssert (source->type == PM_SOURCE_TYPE_UNKNOWN, "source type was not init-ed?"); @@ -561,5 +561,5 @@ if (!source->moments) { source->type = PM_SOURCE_TYPE_STAR; - psAssert (source->mode & noMoments, "why is this source missing moments?"); + psAssert (source->mode & noMoments, "why is this source missing moments?"); Nstar++; continue; @@ -596,36 +596,38 @@ } - // likely defect (too small to be stellar) (push out to 3 sigma) - // low S/N objects which are small are probably stellar - // only set candidate defects if - // XXX these limits are quite arbitrary - if ((sigX < 0.05) || (sigY < 0.05)) { - source->type = PM_SOURCE_TYPE_DEFECT; - source->mode |= PM_SOURCE_MODE_DEFECT; - Ncr ++; - continue; - } - - // likely unsaturated extended source (too large to be stellar) - if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - source->type = PM_SOURCE_TYPE_EXTENDED; - Next ++; - continue; - } - - // the rest are probable stellar objects - starsn_moments->data.F32[starsn_moments->n] = source->moments->SN; - starsn_moments->n ++; - starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN; - starsn_peaks->n ++; - Nstar ++; - - // PSF star (within 1.5 sigma of clump center, S/N > limit) - psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); - if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) { - source->type = PM_SOURCE_TYPE_STAR; - source->mode |= PM_SOURCE_MODE_PSFSTAR; - Npsf ++; - continue; + // The following determinations require the use of moments + if (!(source->mode & noMoments)) { + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // XXX these limits are quite arbitrary + if (sigX < 0.05 || sigY < 0.05) { + source->type = PM_SOURCE_TYPE_DEFECT; + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // likely unsaturated extended source (too large to be stellar) + if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) { + source->type = PM_SOURCE_TYPE_EXTENDED; + Next ++; + continue; + } + + // the rest are probable stellar objects + starsn_moments->data.F32[starsn_moments->n] = source->moments->SN; + starsn_moments->n ++; + starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN; + starsn_peaks->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center, S/N > limit) + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) { + source->type = PM_SOURCE_TYPE_STAR; + source->mode |= PM_SOURCE_MODE_PSFSTAR; + Npsf ++; + continue; + } } @@ -642,8 +644,8 @@ if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); - psFree (stats); - psFree (starsn_peaks); - return false; + psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); + psFree (stats); + psFree (starsn_peaks); + return false; } psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats->min, stats->max); @@ -656,8 +658,8 @@ stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); - psFree (stats); - psFree (starsn_peaks); - return false; + psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats"); + psFree (stats); + psFree (starsn_peaks); + return false; } psLogMsg ("psModules.objects", 3, "SN range (peaks) : %f - %f (%ld)\n", Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.c (revision 24557) @@ -496,4 +496,7 @@ status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); + } if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { @@ -506,4 +509,7 @@ status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); + } } if (xfitname) { @@ -517,4 +523,7 @@ status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); } + if (!strcmp (exttype, "PS1_V2")) { + status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); + } } if (!status) { @@ -944,4 +953,7 @@ sources = pmSourcesRead_CMF_PS1_V1 (file->fits, hdu->header); } + if (!strcmp (exttype, "PS1_V2")) { + sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); + } } Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.h (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.h (revision 24557) @@ -39,4 +39,8 @@ bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname); + bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); @@ -48,4 +52,5 @@ psArray *pmSourcesRead_PS1_CAL_0 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); bool pmSourcesWritePSFs (psArray *sources, char *filename); Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 24557) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 24557) @@ -0,0 +1,625 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.h" + +// panstarrs-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, + psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + psF32 errMag, chisq, apRadius; + psS32 nPix, nDOF; + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + bool status1 = false; + bool status2 = false; + float magOffset = NAN; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + + // if the sequence is defined, write these in seq order; otherwise + // write them in S/N order: + if (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state + // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // If source->seq is -1, source was generated in this analysis. If source->seq is + // not -1, source was read from elsewhere: in the latter case, preserve the source + // ID. source.seq is used instead of source.id since the latter is a const + // generated on Alloc, and would thus be wrong for read in sources. + if (source->seq == -1) { + source->seq = i; + } + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + } else { + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + } + chisq = model->chisq; + nDOF = model->nDOF; + nPix = model->nPix; + apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + chisq = NAN; + nDOF = 0; + nPix = 0; + apRadius = NAN; + errMag = NAN; + } + + float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); + + // distinguish moments measure from window vs S/N > XX ?? + float mxx = source->moments ? source->moments->Mxx : NAN; + float mxy = source->moments ? source->moments->Mxy : NAN; + float myy = source->moments ? source->moments->Myy : NAN; + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + // XXX not sure how to get this : need to load Nimages with weight? + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + + psArrayAdd (table, 100, row); + psFree (row); + + // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not + // subtracted + + // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image + // edge; 3) any pixels in the 3x3 peak region are masked; + } + + psMetadata *header = psMetadataCopy(NULL, tableHeader); + pmSourceMasksHeader(header); + + if (table->n == 0) { + psFitsWriteBlank(fits, header, extname); + psFree(table); + psFree(header); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + assert (modelType > -1); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entriesa + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN; + + pmPSF_AxesToModel (PAR, axes); + + float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + source->peak->flux = peakFlux; + source->peak->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + source->moments = pmMomentsAlloc (); + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + + bool status; + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + assert (radialBinsLower->n == radialBinsUpper->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + // skip source if it is not a ext sourc + // XXX we have two places that extended source parameters are measured: + // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, + // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) + // should we require both? + + pmSource *source = sources->data[i]; + + // skip sources without measurements + if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; + + // XXX I need to split the extended models from the extended aperture measurements + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + if (doPetrosian) { + pmSourcePetrosianValues *petrosian = source->extpars->petrosian; + if (petrosian) { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + } + } + + // Kron measurements + if (doKron) { + pmSourceKronValues *kron = source->extpars->kron; + if (kron) { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); + } + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { + pmSourceIsophotalValues *isophot = source->extpars->isophot; + if (isophot) { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); + } + } + + // Flux Annuli + if (doAnnuli) { + pmSourceAnnuli *annuli = source->extpars->annuli; + if (annuli) { + psVector *fluxVal = annuli->flux; + psVector *fluxErr = annuli->fluxErr; + psVector *fluxVar = annuli->fluxVar; + + for (int j = 0; j < fluxVal->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); + } + } else { + for (int j = 0; j < radialBinsLower->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); + } + } + } + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + + return true; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 24557) @@ -49,5 +49,4 @@ pmCell *cell = NULL; pmReadout *readout = NULL; - pmFPAview *view = pmFPAviewAlloc (0); // first, check if there are any matches to be written @@ -57,6 +56,6 @@ psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSASTRO"); if (!status) { - psError(PS_ERR_UNKNOWN, true, "missing recipe PSASTRO in config data"); - return false; + psError(PS_ERR_UNKNOWN, true, "missing recipe PSASTRO in config data"); + return false; } @@ -65,4 +64,5 @@ psArray *table = psArrayAllocEmpty (0x1000); + pmFPAview *view = pmFPAviewAlloc (0); // this loop selects the matched stars for all chips @@ -70,58 +70,58 @@ psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) continue; - - char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) continue; - // process each of the readouts - // XXX there can only be one readout per chip, right? - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (! readout->data_exists) continue; + // process each of the readouts + // XXX there can only be one readout per chip, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) continue; - // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); - if (rawstars == NULL) continue; + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) continue; - // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) continue; - psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) continue; + psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); - psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); - if (matches == NULL) continue; + psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (matches == NULL) continue; - for (int i = 0; i < matches->n; i++) { - pmAstromMatch *match = matches->data[i]; + for (int i = 0; i < matches->n; i++) { + pmAstromMatch *match = matches->data[i]; - pmAstromObj *raw = rawstars->data[match->raw]; - pmAstromObj *ref = refstars->data[match->ref]; + pmAstromObj *raw = rawstars->data[match->raw]; + pmAstromObj *ref = refstars->data[match->ref]; - psMetadata *row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "RA", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); - psMetadataAdd (row, PS_LIST_TAIL, "DEC", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); - psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); - psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color", ref->Color); - psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); - // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined. + psMetadata *row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "RA", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); + psMetadataAdd (row, PS_LIST_TAIL, "DEC", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color", ref->Color); + psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); + // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined. - psArrayAdd (table, 100, row); - psFree (row); - } - } - } + psArrayAdd (table, 100, row); + psFree (row); + } + } + } } psFree (view); if (table->n == 0) { - psFree(table); - return true; + psFree(table); + return true; } Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceMatch.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceMatch.c (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceMatch.c (revision 24557) @@ -45,4 +45,5 @@ psVector **x, psVector **y, // Coordinate vectors to return psVector **mag, psVector **magErr, // Magnitude and error vectors to return + psVector **indices, // Indices for sources const psArray *sources // Input sources ) @@ -58,13 +59,13 @@ *mag = psVectorRecycle(*mag, numSources, PS_TYPE_F32); *magErr = psVectorRecycle(*magErr, numSources, PS_TYPE_F32); + *indices = psVectorRecycle(*indices, numSources, PS_TYPE_S32); float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources long num = 0; // Number of valid sources for (long i = 0; i < numSources; i++) { pmSource *source = sources->data[i]; // Source of interest - if (!source) continue; - if (source->mode & SOURCE_MASK) continue; - if (!isfinite(source->psfMag)) continue; - if (!isfinite(source->errMag)) continue; - if (source->psfMag > SOURCE_FAINTEST) continue; + if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || + !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) { + continue; + } float xSrc, ySrc; // Coordinates of source @@ -78,5 +79,6 @@ (*y)->data.F32[num] = ySrc; (*mag)->data.F32[num] = source->psfMag; - (*magErr)->data.F32[num] = source->errMag ; + (*magErr)->data.F32[num] = source->errMag; + (*indices)->data.S32[num] = i; num++; } @@ -85,4 +87,5 @@ (*mag)->n = num; (*magErr)->n = num; + (*indices)->n = num; if (*bounds) { @@ -175,6 +178,7 @@ psVector *xImage = NULL, *yImage = NULL; // Coordinates of sources psVector *magImage = NULL, *magErrImage = NULL; // Magnitude and mag - - int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, + psVector *indices = NULL; // Indices for sources + + int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, &indices, sources); // Number of sources @@ -187,5 +191,6 @@ for (int j = 0; j < numSources; j++) { pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); matches->data[j] = match; } @@ -206,5 +211,6 @@ for (int j = 0, k = numMaster; j < numSources; j++, k++) { pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); matches->data[k] = match; } @@ -231,10 +237,12 @@ // Record the match pmSourceMatch *match = matches->data[index]; // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); numMatch++; } else { // Add to the master list pmSourceMatch *match = pmSourceMatchAlloc(); // Match data - pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], + i, indices->data.S32[j]); xMaster->data.F32[numMaster] = xImage->data.F32[j]; yMaster->data.F32[numMaster] = yImage->data.F32[j];