Index: trunk/psModules/src/detrend/pmFringeStats.c =================================================================== --- trunk/psModules/src/detrend/pmFringeStats.c (revision 7836) +++ trunk/psModules/src/detrend/pmFringeStats.c (revision 7875) @@ -3,6 +3,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-07-07 03:26:22 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-07-12 04:44:01 $ * * Copyright 2004 IfA @@ -56,5 +56,5 @@ } -bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image) +bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image, psRandom *random) { PS_ASSERT_PTR_NON_NULL(fringe, false); @@ -65,5 +65,10 @@ // create fringe->nRequested - psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0); + psRandom *rng; + if (random) { + rng = psMemIncrRefCounter(random); + } else { + rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + } fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32); @@ -84,12 +89,142 @@ // generate random points located within image bounds for (int i = 0; i < fringe->nRequested; i++) { - frnd = psRandomUniform(rnd); + frnd = psRandomUniform(rng); xPt[i] = (nX - 2*dX)* frnd + dX; - frnd = psRandomUniform(rnd); + frnd = psRandomUniform(rng); yPt[i] = (nY - 2*dY)* frnd + dY; } + + psFree(rng); + return true; } +bool pmFringeRegionsWriteFits(psFits *fits, psMetadata *header, + const pmFringeRegions *regions, const char *extname) +{ + // Make sure the input is well-behaved + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(regions, false); + psVector *x = regions->x; // The x positions + psVector *y = regions->y; // The y positions + psVector *mask = regions->mask; // The region mask + int numRows = regions->nRequested; // Number of rows in the table + PS_ASSERT_INT_POSITIVE(numRows, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_SIZE(x, numRows, false); + PS_ASSERT_VECTOR_SIZE(y, numRows, false); + if (mask) { + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + PS_ASSERT_VECTOR_SIZE(mask, numRows, false); + } + + // We need to write: + // Scalars: dX, dY, nX, nY + // Vectors: x, y + + psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header + psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width", + regions->dX); + psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height", + regions->dY); + psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x", + regions->nX); + psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y", + regions->nY); + + psArray *table = psArrayAlloc(numRows); // The table + table->n = numRows; + // Translate the vectors into the required format for psFitsWriteTable() + for (long i = 0; i < numRows; i++) { + psMetadata *row = psMetadataAlloc(); + psMetadataAddF32(row, PS_LIST_TAIL, "x", PS_META_REPLACE, "Fringe position in x", x->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "y", PS_META_REPLACE, "Fringe position in y", y->data.F32[i]); + table->data[i] = row; + } + + bool success; // Success of operation + if (!(success = psFitsWriteTable(fits, scalars, table, extname))) { + psError(PS_ERR_IO, false, "Unable to write fringe data to extension %s\n", extname); + } + psFree(scalars); + psFree(table); + + return success; +} + +pmFringeRegions *pmFringeRegionsReadFits(psMetadata *header, const psFits *fits, const char *extname) +{ + PS_ASSERT_PTR_NON_NULL(fits, NULL); + + if (extname && strlen(extname) > 0) { + if (!psFitsMoveExtName(fits, extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname); + return NULL; + } + } else if (!psFitsMoveExtNum(fits, 0, false)) { + psError(PS_ERR_IO, false, "Unable to move to PHU\n"); + return NULL; + } + + psMetadata *headerCopy = psMemIncrRefCounter(header); // Copy of the header, or NULL + + headerCopy = psFitsReadHeader(headerCopy, fits); // The FITS header + if (!headerCopy) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname); + psFree(header); + return NULL; + } + + // Read the scalars from the header + #define READ_SCALAR(SCALAR, NAME) \ + int SCALAR = psMetadataLookupS32(&mdok, header, NAME); \ + if (!mdok || SCALAR <= 0) { \ + psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s.\n", extname); \ + psFree(header); \ + return NULL; \ + } + + // Need to retrieve the scalars: dX, dY, nX, nY + bool mdok = true; // Status of MD lookup + READ_SCALAR(dX, "PSFRNGDX"); + READ_SCALAR(dY, "PSFRNGDY"); + READ_SCALAR(nX, "PSFRNGNX"); + READ_SCALAR(nY, "PSFRNGNY"); + psFree(header); + + // Now the vectors: x, y + psArray *table = psFitsReadTable(fits); // The table + long numRows = table->n; // Number of rows + + pmFringeRegions *regions = pmFringeRegionsAlloc(numRows, dX, dY, nX, nY); // The fringe regions + psVector *x = psVectorAlloc(numRows, PS_TYPE_F32); // x position + psVector *y = psVectorAlloc(numRows, PS_TYPE_F32); // y position + x->n = y->n = numRows; + regions->x = x; + regions->y = y; + + #define READ_REGIONS_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \ + VECTOR->data.TYPE[i] = psMetadataLookup##TYPE(&mdok, row, NAME); \ + if (!mdok) { \ + psError(PS_ERR_IO, true, "Unable to find " #DESCRIPTION " .\n"); \ + psFree(table); \ + psFree(regions); \ + return NULL; \ + } + + // Translate the table into vectors + for (long i = 0; i < numRows; i++) { + psMetadata *row = table->data[i]; // Table row + READ_REGIONS_ROW(x, F32, "x", "x position"); + READ_REGIONS_ROW(y, F32, "y", "y position"); + } + psFree(table); + + return regions; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -129,5 +264,5 @@ if (!fringe->x || !fringe->y) { // create the fringe vectors for this image - pmFringeRegionsCreatePoints(fringe, readout->image); + pmFringeRegionsCreatePoints(fringe, readout->image, NULL); } @@ -167,4 +302,6 @@ psImageStats(median, subImage, subMask, maskVal); sky->data.F32[i][j] = median->sampleMedian; + psFree(subImage); + psFree(subMask); } } @@ -193,9 +330,16 @@ (int)region.y0, (int)region.y1, fPt[i], dfPt[i]); } + psFree(sky); + psFree(median); + psFree(medianSd); return measurements; } -bool pmFringeStatsWriteFits(psFits *fits, const pmFringeStats *fringe, const char *extname) +bool pmFringeStatsWriteFits(psFits *fits, + psMetadata *header, + const pmFringeStats *fringe, + const char *extname + ) { // Make sure the input is well-behaved @@ -204,19 +348,6 @@ pmFringeRegions *regions = fringe->regions; // The fringe regions PS_ASSERT_PTR_NON_NULL(regions, false); - psVector *x = regions->x; // The x positions - psVector *y = regions->y; // The y positions - psVector *mask = regions->mask; // The region mask - int numRows = x->n; // Number of rows in the table - PS_ASSERT_INT_POSITIVE(regions->nRequested, false); - PS_ASSERT_VECTOR_NON_NULL(x, false); - PS_ASSERT_VECTOR_NON_NULL(y, false); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_SIZE(y, numRows, false); - if (mask) { - PS_ASSERT_VECTOR_NON_NULL(mask, false); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); - PS_ASSERT_VECTOR_SIZE(mask, numRows, false); - } + int numRows = regions->nRequested; // Number of rows in the table + PS_ASSERT_INT_POSITIVE(numRows, false); psVector *f = fringe->f; // The fringe measurements psVector *df = fringe->df; // The fringe standard deviatiations @@ -229,17 +360,5 @@ // We need to write: - // Scalars: dX, dY, nX, nY - // Vectors: x, y, mask, f, df - - psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header - psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width", - regions->dX); - psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height", - regions->dY); - psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x", - regions->nX); - psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y", - regions->nY); - + // Vectors: f, df psArray *table = psArrayAlloc(numRows); // The table table->n = numRows; @@ -247,30 +366,27 @@ for (long i = 0; i < numRows; i++) { psMetadata *row = psMetadataAlloc(); - psArraySet(table, i, row); - psMetadataAddF32(row, PS_LIST_TAIL, "x", PS_META_REPLACE, "Fringe position in x", x->data.F32[i]); - psMetadataAddF32(row, PS_LIST_TAIL, "y", PS_META_REPLACE, "Fringe position in y", y->data.F32[i]); - if (mask) { - psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", mask->data.F32[i]); - } else { - psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", 0); - } psMetadataAddF32(row, PS_LIST_TAIL, "f", PS_META_REPLACE, "Fringe measurement", f->data.F32[i]); psMetadataAddF32(row, PS_LIST_TAIL, "df", PS_META_REPLACE, "Fringe stdev", df->data.F32[i]); - psFree(row); // Drop reference - } - - bool success; // Success of operation - if (!(success = psFitsWriteTable(fits, scalars, table, extname))) { + table->data[i] = row; + } + + if (!psFitsWriteTable(fits, header, table, extname)) { psError(PS_ERR_IO, false, "Unable to write fringe data to extension %s\n", extname); - } - psFree(scalars); + psFree(table); + return false; + } + psFree(table); - - return success; -} - -pmFringeStats *pmFringeStatsReadFits(const psFits *fits, const char *extname) + return true; +} + +pmFringeStats *pmFringeStatsReadFits(psMetadata *header, const psFits *fits, const char *extname, + pmFringeRegions *regions) { PS_ASSERT_PTR_NON_NULL(fits, NULL); + PS_ASSERT_PTR_NON_NULL(regions, NULL); + PS_ASSERT_INT_POSITIVE(regions->nRequested, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(regions->x, regions->y, NULL); + PS_ASSERT_VECTOR_SIZE(regions->x, regions->nRequested, NULL); if (extname && strlen(extname) > 0) { @@ -284,44 +400,23 @@ } - psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header - if (!header) { + psMetadata *headerCopy = psMemIncrRefCounter(header); // Copy of the header, or NULL + + headerCopy = psFitsReadHeader(headerCopy, fits); // The FITS header + if (!headerCopy) { psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname); + psFree(headerCopy); return NULL; } - - // Read the scalars from the header - #define READ_SCALAR(SCALAR, NAME) \ - int SCALAR = psMetadataLookupS32(&mdok, header, NAME); \ - if (!mdok || SCALAR <= 0) { \ - psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s.\n", extname); \ - return NULL; \ - } - - // Need to retrieve the scalars: dX, dY, nX, nY - bool mdok = true; // Status of MD lookup - READ_SCALAR(dX, "PSFRNGDX"); - READ_SCALAR(dY, "PSFRNGDX"); - READ_SCALAR(nX, "PSFRNGDX"); - READ_SCALAR(nY, "PSFRNGDX"); - psFree(header); - - // Now the vectors: x, y, mask, f, df + psFree(headerCopy); + + // Now the vectors: f, df psArray *table = psFitsReadTable(fits); // The table long numRows = table->n; // Number of rows - pmFringeRegions *regions = pmFringeRegionsAlloc(numRows, dX, dY, nX, nY); // The fringe regions - psVector *x = psVectorAlloc(numRows, PS_TYPE_F32); // x position - psVector *y = psVectorAlloc(numRows, PS_TYPE_F32); // y position - psVector *mask = psVectorAlloc(numRows, PS_TYPE_U8); // mask - x->n = y->n = mask->n = numRows; - regions->x = x; - regions->y = y; - regions->mask = mask; pmFringeStats *fringes = pmFringeStatsAlloc(regions); // The fringe measurements - psFree(regions); // Drop reference psVector *f = fringes->f; // fringe measurement - psVector *df = fringes->df; // fringe stdev - - #define READ_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \ + psVector *df = fringes->df; // fringe stdev + + #define READ_STATS_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \ VECTOR->data.TYPE[i] = psMetadataLookup##TYPE(&mdok, row, NAME); \ if (!mdok) { \ @@ -333,11 +428,9 @@ // Translate the table into vectors + bool mdok; // Status of MD lookup for (long i = 0; i < numRows; i++) { psMetadata *row = table->data[i]; // Table row - READ_ROW(x, F32, "x", "x position"); - READ_ROW(y, F32, "y", "y position"); - READ_ROW(mask, U8, "mask", "mask value"); - READ_ROW(f, F32, "f", "fringe measurement"); - READ_ROW(df, F32, "df", "fringe standard deviation"); + READ_STATS_ROW(f, F32, "f", "fringe measurement"); + READ_STATS_ROW(df, F32, "df", "fringe standard deviation"); } psFree(table);