Index: /branches/pap_branch_051214/psModules/src/astrom/Makefile.am =================================================================== --- /branches/pap_branch_051214/psModules/src/astrom/Makefile.am (revision 5794) +++ /branches/pap_branch_051214/psModules/src/astrom/Makefile.am (revision 5795) @@ -4,9 +4,25 @@ libpsmoduleastrom_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmoduleastrom_la_SOURCES = \ - pmAstrometry.c \ - pmAstrometryObjects.c + pmAstrometry.c \ + pmAstrometryObjects.c \ + pmChipMosaic.c \ + pmFPAConceptsGet.c \ + pmFPAConceptsSet.c \ + pmFPAConstruct.c \ + pmFPARead.c \ + pmFPAWrite.c \ + pmReadout.c \ + psAdditionals.c psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ - pmAstrometry.h \ - pmAstrometryObjects.h + pmAstrometry.h \ + pmAstrometryObjects.h \ + pmChipMosaic.h \ + pmFPAConceptsGet.h \ + pmFPAConceptsSet.h \ + pmFPAConstruct.h \ + pmFPARead.h \ + pmFPAWrite.h \ + pmReadout.h \ + psAdditionals.h Index: /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.c =================================================================== --- /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.c (revision 5794) +++ /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.c (revision 5795) @@ -10,9 +10,9 @@ * XXX: We should review the extent of the warning messages on these functions * when the transformations are not successful. -* +* * XXX: Should we implement non-linear cell->chip transforms? -* -* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-08 00:00:57 $ +* +* @version $Revision: 1.11.2.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-17 03:18:38 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -60,19 +60,4 @@ psFree(readout->mask); psFree(readout->weight); - // - // XXX: Not sure if this is the right way to do things. Currently the psListAdd() - // increase the memory reference counter to the list data. So, we - // iterate through the list, and decrement the reference counters. - // - if (1) { - if ((readout->bias != NULL) && (readout->bias->head != NULL)) { - psListElem *tmpElem = (psListElem *) readout->bias->head; - while (NULL != tmpElem) { - psMemDecrRefCounter((psImage *) tmpElem->data); - tmpElem = tmpElem->next; - } - } - } - psFree(readout->bias); psFree(readout->analysis); psFree(readout->parent); @@ -88,4 +73,5 @@ psFree(cell->concepts); psFree(cell->analysis); + psFree(cell->camera); // // Set the parent to NULL in all cell->readouts before psFree(cell->readouts) @@ -101,5 +87,5 @@ psFree(cell->readouts); psFree(cell->parent); - psFree(cell->private); + psFree(cell->hdu); } } @@ -125,5 +111,5 @@ psFree(chip->cells); psFree(chip->parent); - psFree(chip->private); + psFree(chip->hdu); } } @@ -151,22 +137,33 @@ } psFree(fpa->chips); - psFree(fpa->private); + psFree(fpa->hdu); psFree(fpa->phu); } } +void p_pmHDUFree(p_pmHDU *hdu) +{ + if (hdu) { + psFree(hdu->extname); + psFree(hdu->header); + psFree(hdu->images); + psFree(hdu->masks); + psFree(hdu->weights); + } +} + // XXX: Verify these default values for row0, col0, rowBins, colBins +// PAP: These values may disappear in the future in favour of values in parent->concepts? pmReadout *pmReadoutAlloc(pmCell *cell) { pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout)); - tmpReadout->col0 = -1; - tmpReadout->row0 = -1; - tmpReadout->colBins = -1; - tmpReadout->rowBins = -1; + tmpReadout->col0 = 0; + tmpReadout->row0 = 0; + tmpReadout->colBins = 0; + tmpReadout->rowBins = 0; tmpReadout->image = NULL; tmpReadout->mask = NULL; tmpReadout->weight = NULL; - tmpReadout->bias = NULL; tmpReadout->analysis = psMetadataAlloc(); tmpReadout->parent = cell; @@ -179,13 +176,14 @@ // XXX: Verify these default values for row0, col0. +// PAP: These values may disappear in the future in favour of values in the "concepts"? pmCell *pmCellAlloc( pmChip *chip, - psMetadata *cameradata, - psString name) + psMetadata *cameraData, + const char *name) { pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell)); - tmpCell->col0 = -1; - tmpCell->row0 = -1; + tmpCell->col0 = 0; + tmpCell->row0 = 0; tmpCell->toChip = NULL; tmpCell->toFPA = NULL; @@ -196,13 +194,13 @@ psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n"); } - tmpCell->camera = cameradata; - tmpCell->analysis = NULL; + tmpCell->camera = psMemIncrRefCounter(cameraData); + tmpCell->analysis = psMetadataAlloc(); tmpCell->readouts = psArrayAlloc(0); - tmpCell->parent = chip; + tmpCell->parent = psMemIncrRefCounter(chip); if (chip != NULL) { chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell); } - tmpCell->valid = false; - tmpCell->private = NULL; + tmpCell->valid = true; // All cells are valid by default + tmpCell->hdu = NULL; psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree); @@ -211,12 +209,13 @@ // XXX: Verify these default values for row0, col0. +// PAP: row0, col0 may disappear in the future in favour of storing values in the "concepts". pmChip *pmChipAlloc( pmFPA *fpa, - psString name) + const char *name) { pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip)); - tmpChip->col0 = -1; - tmpChip->row0 = -1; + tmpChip->col0 = 0; + tmpChip->row0 = 0; tmpChip->toFPA = NULL; tmpChip->fromFPA = NULL; @@ -226,12 +225,12 @@ psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME to metadata.\n"); } - tmpChip->analysis = NULL; + tmpChip->analysis = psMetadataAlloc(); tmpChip->cells = psArrayAlloc(0); - tmpChip->parent = fpa; + tmpChip->parent = psMemIncrRefCounter(fpa); if (fpa != NULL) { fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip); } - tmpChip->valid = false; - tmpChip->private = NULL; + tmpChip->valid = true; // Work on all chips, by default + tmpChip->hdu = NULL; psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree); @@ -248,11 +247,25 @@ tmpFPA->concepts = psMetadataAlloc(); tmpFPA->analysis = NULL; - tmpFPA->camera = camera; + tmpFPA->camera = psMemIncrRefCounter((psPtr)camera); tmpFPA->chips = psArrayAlloc(0); - tmpFPA->private = NULL; + tmpFPA->hdu = NULL; tmpFPA->phu = NULL; psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree); return(tmpFPA); +} + +p_pmHDU *p_pmHDUAlloc(const char *extname) +{ + p_pmHDU *hdu = psAlloc(sizeof(p_pmHDU)); + psMemSetDeallocator(hdu, (psFreeFunc)p_pmHDUFree); + + hdu->extname = psStringCopy(extname); + hdu->header = NULL; + hdu->images = NULL; + hdu->masks = NULL; + hdu->weights = NULL; + + return hdu; } @@ -800,4 +813,29 @@ } + +/***************************************************************************** + *****************************************************************************/ + +// Set cells within a chip to a specified validity +static bool setCellsValid(const pmChip *chip, // Chip of interest + bool valid // Valid? + ) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + + psArray *cells = chip->cells; // Component cells + if (! cells) { + return false; + } + for (int i = 0; i < cells->n; i++) { + pmCell *tmpCell = cells->data[i]; // Cell of interest + if (tmpCell) { + tmpCell->valid = valid; + } + } + + return true; +} + /***************************************************************************** *****************************************************************************/ @@ -807,23 +845,26 @@ { PS_ASSERT_PTR_NON_NULL(fpa, false); - if ((fpa->chips == NULL) || (chipNum >= fpa->chips->n)) { + + psArray *chips = fpa->chips; // Component chips + if ((chips == NULL) || (chipNum >= chips->n)) { return(false); } - psBool rc = true; - - for (psS32 i = 0 ; i < fpa->chips->n ; i++) { - pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; + + for (int i = 0 ; i < chips->n ; i++) { + pmChip *tmpChip = (pmChip *) chips->data[i]; if (tmpChip == NULL) { - rc = false; + continue; + } + if (i == chipNum) { + tmpChip->valid = true; + setCellsValid(tmpChip, true); } else { - if (i == chipNum) { - tmpChip->valid = true; - } else { - tmpChip->valid = false; - } - } - } - - return(rc); + tmpChip->valid = false; + setCellsValid(tmpChip, false); + } + + } + + return true; } @@ -831,9 +872,9 @@ /***************************************************************************** XXX: The SDRS is ambiguous on a few things: - Whether or not the other chips should be set valid=true. - Should we return the number of chip valid=true before or after they're set, + Whether or not the other chips should be set valid=true. [PAP: No] + Should we return the number of chip valid=true before or after they're set, [PAP: After] *****************************************************************************/ /** - * + * * pmFPAExcludeChip shall set valid to false only for the specified chip * number (chipNum). In the event that the specified chip number does not exist @@ -841,5 +882,5 @@ * The function shall return the number of chips within the fpa that have valid * set to true. - * + * */ int pmFPAExcludeChip( @@ -849,21 +890,22 @@ PS_ASSERT_PTR_NON_NULL(fpa, false); - if (fpa->chips == NULL) { + psArray *chips = fpa->chips; // Component chips + if (chips == NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n"); return(0); } - if ((chipNum >= fpa->chips->n) || (NULL == (pmChip *) fpa->chips->data[chipNum])) { + if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum); return(0); } - psS32 numChips = 0; - for (psS32 i = 0 ; i < fpa->chips->n ; i++) { - pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; + int numChips = 0; // Number of valid chips + for (int i = 0 ; i < chips->n ; i++) { + pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest if (tmpChip != NULL) { if (i == chipNum) { tmpChip->valid = false; - } else { - tmpChip->valid = true; + setCellsValid(tmpChip, false); // Wipe out the cell as well + } else if (tmpChip->valid) { numChips++; } Index: /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.h =================================================================== --- /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/astrom/pmAstrometry.h (revision 5795) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-08 00:00:57 $ +* @version $Revision: 1.5.2.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-17 03:18:38 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -30,11 +30,13 @@ { const char *extname; // Extension name, if it corresponds to this level - psArray *pixels; // The pixel data, if it corresponds to this level psMetadata *header; // The FITS header, if it corresponds to this level + psArray *images; // The pixel data, if it corresponds to this level + psArray *masks; // The mask data, if it corresponds to this level + psArray *weights; // The weight data, if it corresponds to this level } p_pmHDU; /** Focal plane data structure - * + * * A focal plane consists of one or more chips (according to the number of * pieces of contiguous silicon). It contains metadata containers for the @@ -49,5 +51,5 @@ * transformation which will be derived from numerically inverting the forward * transformation. - * + * */ typedef struct @@ -62,5 +64,5 @@ const psMetadata *camera; ///< Camera configuration psArray *chips; ///< The chips - p_pmHDU *private; ///< FITS data + p_pmHDU *hdu; ///< FITS data psMetadata *phu; ///< Primary Header } @@ -68,5 +70,5 @@ /** Chip data structure - * + * * A chip consists of one or more cells (according to the number of amplifiers * on the device). The chip contains metadata containers for the concepts and @@ -79,5 +81,5 @@ * inverting the forward transformation. A boolean indicates whether the chip is * of interest, allowing it to be excluded from analysis. - * + * */ typedef struct @@ -95,5 +97,5 @@ pmFPA *parent; ///< Parent FPA bool valid; ///< Do we bother about reading and working with this chip? - p_pmHDU *private; ///< FITS data + p_pmHDU *hdu; ///< FITS data } pmChip; @@ -124,5 +126,5 @@ pmChip *parent; ///< Parent chip bool valid; ///< Do we bother about reading and working with this cell? - p_pmHDU *private; ///< FITS data + p_pmHDU *hdu; ///< FITS data } pmCell; @@ -141,4 +143,5 @@ { // Position on the cell + // XXX: These may be removed in the future; use parent->concepts instead? int col0; ///< Offset from the left of chip. int row0; ///< Offset from the bottom of chip. @@ -149,5 +152,4 @@ psImage *mask; ///< Mask of input image psImage *weight; ///< Weight of input image - psList *bias; ///< List of bias section (sub-)images psMetadata *analysis; ///< Readout-level analysis metadata pmCell *parent; ///< Parent cell @@ -180,11 +182,11 @@ */ pmCell *pmCellAlloc( - pmChip *chip, ///< Parent chip - psMetadata *cameradata, - psString name + pmChip *chip, ///< Parent chip + psMetadata *cameraData, ///< Camera data + const char *name ///< Name of cell ); /** Allocates a pmChip - * + * * The constructor shall make an empty pmChip. If the parent fpa is not NULL, * the parent link is made and the chip shall be placed in the parent's array @@ -196,16 +198,15 @@ */ pmChip *pmChipAlloc( - pmFPA *fpa, - psString name - + pmFPA *fpa, ///< FPA to which the chip belongs + const char *name ///< Name of chip ); /** Allocates a pmFPA - * + * * The constructor shall make an empty pmFPA. The chips array shall be * allocated with a zero size, the camera and db pointers set to the values * provided, and the concepts metadata constructed. All other pointers in the * structure shall be initialized to NULL. - * + * */ pmFPA *pmFPAAlloc( @@ -213,12 +214,21 @@ ); +/** Allocates a p_pmHDU + * + * XXX: More detailed description + * + * @return p_pmHDU* newly allocated p_pmHDU + */ +p_pmHDU *p_pmHDUAlloc(const char *extname // Extension name + ); + /** Verify parent links. - * + * * This function checks the validity of the parent links in the FPA hierarchy. * If a parent link is not set (or not set correctly), it is corrected, and the * function shall return false. If all the parent pointers were correct, the * function shall return true. - * + * */ bool pmFPACheckParents( @@ -229,10 +239,38 @@ /** FUNC DESC - * - * - * - * - */ - + * + * + * + * + */ + + + +/** + * + * pmFPASelectChip shall set valid to true for the specified chip number + * (chipNum), and all other chips shall have valid set to false. In the event + * that the specified chip number does not exist within the fpa, the function + * shall return false. + * + */ +bool pmFPASelectChip( + pmFPA *fpa, + int chipNum +); + +/** + * + * pmFPAExcludeChip shall set valid to false only for the specified chip + * number (chipNum). In the event that the specified chip number does not exist + * within the fpa, the function shall generate a warning, and perform no action. + * The function shall return the number of chips within the fpa that have valid + * set to true. + * + */ +int pmFPAExcludeChip( + pmFPA *fpa, + int chipNum +); @@ -423,31 +461,4 @@ ); -/** - * - * pmFPASelectChip shall set valid to true for the specified chip number - * (chipNum), and all other chips shall have valid set to false. In the event - * that the specified chip number does not exist within the fpa, the function - * shall return false. - * - */ -bool pmFPASelectChip( - pmFPA *fpa, - int chipNum -); - -/** - * - * pmFPAExcludeChip shall set valid to false only for the specified chip - * number (chipNum). In the event that the specified chip number does not exist - * within the fpa, the function shall generate a warning, and perform no action. - * The function shall return the number of chips within the fpa that have valid - * set to true. - * - */ -int pmFPAExcludeChip( - pmFPA *fpa, - int chipNum -); - #endif // #ifndef PS_ASTROMETRY_H Index: /branches/pap_branch_051214/psModules/src/camera/Makefile.am =================================================================== --- /branches/pap_branch_051214/psModules/src/camera/Makefile.am (revision 5794) +++ /branches/pap_branch_051214/psModules/src/camera/Makefile.am (revision 5795) @@ -3,7 +3,20 @@ libpsmodulecamera_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) libpsmodulecamera_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmodulecamera_la_SOURCES = +libpsmodulecamera_la_SOURCES = \ + pmChipMosaic.c \ + pmFPAConceptsGet.c \ + pmFPAConceptsSet.c \ + pmFPAConstruct.c \ + pmFPAMorph.c \ + pmFPARead.c \ + pmFPAWrite.c psmoduleincludedir = $(includedir) -psmoduleinclude_HEADERS = - +psmoduleinclude_HEADERS = \ + pmChipMosaic.h \ + pmFPAConceptsGet.h \ + pmFPAConceptsSet.h \ + pmFPAConstruct.h \ + pmFPAMorph.h \ + pmFPARead.h \ + pmFPAWrite.h Index: /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.c =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.c (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.c (revision 5795) @@ -1,2 +1,8 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + /** @file pmFlatField.c * @@ -18,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-18 19:43:14 $ + * @version $Revision: 1.4.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -30,5 +36,5 @@ #include #include -#include +#include #include "pslib.h" @@ -36,32 +42,10 @@ #include "pmMaskBadPixels.h" #include "pmFlatFieldErrors.h" -#include "pmSubtractBias.h" -// XXX: This should be removed when the autoconf stuff handles psConstants.h correctly. -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ -bool pmFlatField( - pmReadout *in, - const pmReadout *flat) +bool pmFlatField(pmReadout *in, pmReadout *mask, const pmReadout *flat) { - // XXX: Use the proper image and readout asserts. - PS_ASSERT_PTR_NON_NULL(in, true); - PS_ASSERT_PTR_NON_NULL(in->image, false); - PS_ASSERT_PTR_NON_NULL(in->mask, false); - PS_ASSERT_PTR_NON_NULL(flat, false); - PS_ASSERT_PTR_NON_NULL(flat->image, false); - if (in == NULL) - printf("XXX: NULL\n"); - // XXX: Not sure if this is correct. Must consult with IfA. - PS_ASSERT_PTR_NON_NULL(in->mask, false); - - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } + PS_ASSERT_PTR_NON_NULL(mask, false); int i = 0; int j = 0; @@ -71,25 +55,37 @@ psElemType flatType; psElemType maskType; + psImage *inImage = NULL; psImage *inMask = NULL; psImage *flatImage = NULL; - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); + + // Check for nulls + if (in == NULL) { + return true; // Readout may not have data in it + } else if(flat==NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT); + return false; + } + + inImage = in->image; flatImage = flat->image; - inMask = in->mask; + if (inImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE); + return false; + } else if(flatImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE); + return false; + } + inMask = mask->image; // Check input image and its mask are not larger than flat image - if (trimmedImg == NULL) - printf("XXX: 00\n"); - if (flatImage == NULL) - printf("XXX 01\n"); - - if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) { + if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE, - trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols); + inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols); return false; } @@ -102,6 +98,6 @@ // Determine total offset based on image offset with chip offset - totOffCol = trimmedImg->col0 + in->col0; - totOffRow = trimmedImg->row0 + in->row0; + totOffCol = inImage->col0 + in->col0; + totOffRow = inImage->row0 + in->row0; // Check that offsets are within image limits @@ -111,8 +107,8 @@ totOffRow, totOffCol, flatImage->numRows, flatImage->numCols); return false; - } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) { + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols); + totOffRow, totOffCol, inImage->numRows, inImage->numCols); return false; } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { @@ -124,5 +120,5 @@ // Check for incorrect types - inType = trimmedImg->type.type; + inType = inImage->type.type; flatType = flatImage->type.type; maskType = inMask->type.type; @@ -153,6 +149,6 @@ case PS_TYPE_##TYPE: \ /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */ \ - for(j = totOffRow; j < trimmedImg->numRows; j++) { \ - for(i = totOffCol; i < trimmedImg->numCols; i++) { \ + for(j = totOffRow; j < inImage->numRows; j++) { \ + for(i = totOffCol; i < inImage->numCols; i++) { \ if(flatImage->data.TYPE[j][i] <= 0.0) { \ /* Negative or zero flat pixels shall be masked in input image as PM_MASK_FLAT */ \ @@ -162,9 +158,9 @@ } \ } \ - for(j = totOffRow; j < trimmedImg->numRows; j++) { \ - for(i = totOffCol; i < trimmedImg->numCols; i++) { \ + for(j = totOffRow; j < inImage->numRows; j++) { \ + for(i = totOffCol; i < inImage->numCols; i++) { \ if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) { \ /* Module shall divide the input image by the flat-fielded image */ \ - trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ + inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ } \ } \ Index: /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.h =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmFlatField.h (revision 5795) @@ -1,2 +1,8 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + /** @file pmFlatField.h * @@ -18,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-20 23:06:24 $ + * @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -37,4 +43,5 @@ bool pmFlatField( pmReadout *in, ///< Readout with input image + pmReadout *mask, ///< Input image mask const pmReadout *flat ///< Readout with flat image ); Index: /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.c (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.c (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixels.c * @@ -19,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-18 19:43:14 $ + * @version $Revision: 1.3.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -31,157 +36,77 @@ #include #include -#include +#include #include "pmMaskBadPixels.h" #include "pmMaskBadPixelsErrors.h" -#include "pmSubtractBias.h" -//XXX: REmove, autoconf is broken. -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ +bool pmMaskBadPixels(pmReadout *in, const psImage *mask, unsigned int maskVal, float sat, + unsigned int growVal, int grow) +{ + int i = 0; + int j = 0; + int jj = 0; + int ii = 0; + int rRound = 0; + int rowMin = 0; + int rowMax = 0; + int colMin = 0; + int colMax = 0; + int totOffCol = 0; + int totOffRow = 0; + float r = 0.0f; + psElemType inType; + psElemType maskType; + psImage *inImage = NULL; + psImage *inMask = NULL; -/****************************************************************************** -GrowPixel(inMask, row, col, radius, maskVal): This private routine takes an -input image mask and a pixel location, then sets (logical or) all pixels with -parameter radius if that pixel to maskVal parameter. -*****************************************************************************/ -psBool GrowPixel( - psImage *inMask, - psS32 col, - psS32 row, - psS32 radius, - psU32 maskVal) -{ - psS32 rowMin = PS_MAX(row-radius, 0); - psS32 rowMax = PS_MIN(row+radius+1, inMask->numRows); - psS32 colMin = PS_MAX(col-radius, 0); - psS32 colMax = PS_MIN(col+radius+1, inMask->numCols); - psF32 squareRadius = PS_SQR(radius); + // Check for nulls + if (in == NULL) { + return true; // Readout may not have data in it + } else if(mask==NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE); + return false; + } - - for(psS32 i=rowMin; idata.U8[i][j] |= maskVal; - } - } - } - return(true); -} - -/****************************************************************************** -GetRadius(inImg, col, row, sat, growVal, grow): This private routine takes an -input image and pixel location and determines what radius that pixel must grow -by. - -//XXX: Inline this or macro it for speed. - *****************************************************************************/ -static psS32 GetRadius( - psImage *inImg, - psS32 col, - psS32 row, - psF32 sat, - psU32 growVal, - psS32 grow) -{ - psS32 growRadius = 0; - if (inImg->type.type == PS_TYPE_F32) { - if(inImg->data.F32[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.F32[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else if (inImg->type.type == PS_TYPE_S32) { - if(inImg->data.S32[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.S32[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else if (inImg->type.type == PS_TYPE_U16) { - if(inImg->data.U16[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.U16[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE, - inImg->type.type); - } - return(growRadius); -} - - -bool pmMaskBadPixels( - pmReadout *in, - const pmReadout *mask, - unsigned int maskVal, - float sat, - unsigned int growVal, - int grow) -{ - // XXX: Review this code then put all asserts at the top. - PS_ASSERT_PTR_NON_NULL(in, true); - PS_ASSERT_PTR_NON_NULL(in->image, false); - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - PS_ASSERT_PTR_NON_NULL(mask, false); - int i = 0; - int j = 0; - int totOffCol = 0; - int totOffRow = 0; - psElemType inType; - psElemType maskType; - psImage *inMask = NULL; - - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - if(in->mask == NULL) { - in->mask = psImageAlloc(trimmedImg->numCols, trimmedImg->numRows, PS_TYPE_MASK); - PS_IMAGE_SET_U8(in->mask, 0); + inImage = in->image; + if (inImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE); + return false; + } else if(in->mask == NULL) { + in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK); + memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK)); } inMask = in->mask; // Check input image and its mask are not larger than mask - if(trimmedImg->numRows > mask->image->numRows || trimmedImg->numCols > mask->image->numCols) { + if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE, - trimmedImg->numRows, trimmedImg->numCols, mask->image->numRows, mask->image->numCols); + inImage->numRows, inImage->numCols, mask->numRows, mask->numCols); return false; - } else if(inMask->numRows > mask->image->numRows || inMask->numCols > mask->image->numCols) { + } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE, - inMask->numRows, inMask->numCols, mask->image->numRows, mask->image->numCols); + inMask->numRows, inMask->numCols, mask->numRows, mask->numCols); return false; } // Determine total offset based on image offset with chip offset - totOffCol = trimmedImg->col0 + in->col0; - totOffRow = trimmedImg->row0 + in->row0; + totOffCol = inImage->col0 + in->col0; + totOffRow = inImage->row0 + in->row0; // Check that offsets are within image limits - if(totOffRow>=mask->image->numRows || totOffCol>=mask->image->numCols) { + if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE, - totOffRow, totOffCol, mask->image->numRows, mask->image->numCols); + totOffRow, totOffCol, mask->numRows, mask->numCols); return false; - } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) { + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols); + totOffRow, totOffCol, inImage->numRows, inImage->numCols); return false; } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { @@ -193,6 +118,6 @@ // Check for incorrect types - inType = trimmedImg->type.type; - maskType = mask->image->type.type; + inType = inImage->type.type; + maskType = mask->type.type; if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { psError( PS_ERR_BAD_PARAMETER_TYPE, true, @@ -207,34 +132,57 @@ } - // - // This is the main loop which examines each pixel and masks pixels if - // A: The mask matches the maskValue - // B: The pixel is larger than the saturation value - // C: The pixel equals the grow value (in which case a circle is masked) - // - for(i=totOffRow; inumRows; i++) { - for(j=totOffCol; jnumCols; j++) { - // - // A: Pixels which satisfy maskVal shall be masked - // - if (mask->image->data.U8[i][j] == maskVal) { - in->mask->data.U8[i][j] |= maskVal; - } + // Macro for all PS types + #define PM_BAD_PIXELS(TYPE) \ +case PS_TYPE_##TYPE: \ + for(j=totOffRow; jnumRows; j++) { \ + for(i=totOffCol; inumCols; i++) { \ + \ + /* Pixels with flux greater than sat shall be masked */ \ + if(inImage->data.TYPE[j][i] > sat) { \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT; \ + } \ + \ + /* Pixels which satisfy maskVal shall be masked */ \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal); \ + \ + /* Pixels which satisfy growVal and within the grow radius shall be masked */ \ + if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) { \ + rowMin = MAX(j-grow, 0); \ + rowMax = MIN(j+grow+1, inImage->numRows); \ + colMin = MAX(i-grow, 0); \ + colMax = MIN(i+grow+1, inImage->numCols); \ + for(jj=rowMin; jjdata.PS_TYPE_MASK_DATA[jj][ii] |= \ + (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal); \ + } \ + } \ + } \ + } \ + } \ + } \ + break; - // - // We first determine how much we need to grow by and store this in - // growRadius. - // - psS32 growRadius = GetRadius(trimmedImg, j, i, sat, growVal, grow); - - // - // Grow the pixel mask if necessary. - // - if (growRadius != 0) { - GrowPixel(in->mask, j, i, growRadius, maskVal); - } - } + // Switch to call bad pixel masking macro defined above + switch(inType) { + PM_BAD_PIXELS(U8); + PM_BAD_PIXELS(U16); + PM_BAD_PIXELS(U32); + PM_BAD_PIXELS(U64); + PM_BAD_PIXELS(S8); + PM_BAD_PIXELS(S16); + PM_BAD_PIXELS(S32); + PM_BAD_PIXELS(S64); + PM_BAD_PIXELS(F32); + PM_BAD_PIXELS(F64); + default: + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED, + inType); } - return true; + return false; } Index: /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixels.h (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixels.h * @@ -19,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-15 20:09:03 $ + * @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -36,5 +41,4 @@ } pmMaskValue; -// XXX: Use PS_MIN, PS_MAX /** Macro to find maximum of two numbers */ #define MAX(A,B)((A)>=(B)?(A):(B)) @@ -54,5 +58,5 @@ bool pmMaskBadPixels( pmReadout *in, ///< Readout containing input image data. - const pmReadout *mask, ///< Mask data to be added to readout mask data. + const psImage *mask, ///< Mask data to be added to readout mask data. unsigned int maskVal, ///< Mask value to determine what to add to input mask. float sat, ///< Saturation limit to mask bad pixels. Index: /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixelsErrors.h =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixelsErrors.h * @@ -7,6 +12,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-28 20:43:52 $ + * @version $Revision: 1.1.16.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -23,5 +28,5 @@ * $2 The error text (rest of the line in pmMaskBadPixelsErrors.h) * $n The order of the source line in pmMaskBadPixelsErrors.h (comments excluded) - * + * * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. */ Index: /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.c =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.c (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.c (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmNonLinear.c * @@ -5,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-05 20:49:40 $ + * @version $Revision: 1.5.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,22 +27,14 @@ #include "pmNonLinear.h" -#include "pmSubtractBias.h" -// XXX: Remove, autoconf must be -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ /****************************************************************************** pmNonLinearityLookup(): This routine will take an pmReadout image as input and a 1-D polynomial. For each pixel in the input image, the polynomial will -be evaluated at that pixels value, and the image pixel will then be set -to +be evaluated at that pixels value, and the image pixel will then be set to that value. -*****************************************************************************/ + *****************************************************************************/ -pmReadout *pmNonLinearityPolynomial( - pmReadout *inputReadout, - const psPolynomial1D *input1DPoly) +pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, + const psPolynomial1D *input1DPoly) { PS_ASSERT_PTR_NON_NULL(inputReadout, NULL); @@ -45,18 +42,11 @@ PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL); - PS_WARN_PTR_NON_NULL(inputReadout->parent); - if (inputReadout->parent != NULL) { - PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); + psS32 i; + psS32 j; - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly, - trimmedImg->data.F32[i][j]); + for (i=0;iimage->numRows;i++) { + for (j=0;jimage->numCols;j++) { + inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]); } } @@ -71,48 +61,68 @@ inFluxe, and the corresponding value in outFlux. The image pixel will then be set to the value from outFlux. - -XXX: Must assert that filename exists. This should probably happen in -the lookup files. *****************************************************************************/ -pmReadout *pmNonLinearityLookup( - pmReadout *inputReadout, - const char *filename -) +pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, + const psVector *inFlux, + const psVector *outFlux) { PS_ASSERT_PTR_NON_NULL(inputReadout,NULL); PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL); PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); - PS_WARN_PTR_NON_NULL(inputReadout->parent); - if (inputReadout->parent != NULL) { - PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts); + PS_ASSERT_PTR_NON_NULL(inFlux,NULL); + if (inFlux->n < 2) { + psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements. Returning inputReadout image."); + return(inputReadout); } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); + PS_ASSERT_PTR_NON_NULL(outFlux,NULL); + psS32 tableSize = inFlux->n; + if (inFlux->n != outFlux->n) { + tableSize = PS_MIN(inFlux->n, outFlux->n); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n); + } + PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL); - psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0); - psS32 numLines = psLookupTableRead(tmpLT); + psS32 i; + psS32 j; + psS32 binNum; + psScalar x; psS32 numPixels = 0; - if (numLines < 2) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Lookup Table is too small. Returning original pmReadout.\n"); - } else { - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1); - if (!isnan(tmpD)) { - trimmedImg->data.F32[i][j] = tmpD; - } else { - numPixels++; - } + psF32 slope; + + x.type.type = PS_TYPE_F32; + for (i=0;iimage->numRows;i++) { + for (j=0;jimage->numCols;j++) { + x.data.F32 = inputReadout->image->data.F32[i][j]; + binNum = p_psVectorBinDisect((psVector *)inFlux, &x); + + if (binNum == -2) { + // We get here if x is below the table lookup range. + inputReadout->image->data.F32[i][j] = outFlux->data.F32[0]; + numPixels++; + + } else if (binNum == -1) { + // We get here if x is above the table lookup range. + inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1]; + numPixels++; + + } else if (binNum < -2) { + // We get here if there was some other problem. + psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect(). Returning inputReadout image."); + return(inputReadout); + numPixels++; + } else { + // Perform linear interpolation. + slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) / + (inFlux->data.F32[binNum+1] - inFlux->data.F32[binNum]); + inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] + + ((x.data.F32 - inFlux->data.F32[binNum]) * slope); } } - if (numPixels > 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); - } } - + if (numPixels > 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); + } return(inputReadout); } Index: /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.h =================================================================== --- /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/detrend/pmNonLinear.h (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmNonLinear.h * @@ -5,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-20 23:06:24 $ + * @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,13 +23,10 @@ #include "pmAstrometry.h" -pmReadout *pmNonLinearityPolynomial( - pmReadout *in, - const psPolynomial1D *coeff -); +pmReadout *pmNonLinearityPolynomial(pmReadout *in, + const psPolynomial1D *coeff); -pmReadout *pmNonLinearityLookup( - pmReadout *in, - const char *filename -); +pmReadout *pmNonLinearityLookup(pmReadout *in, + const psVector *inFlux, + const psVector *outFlux); #endif Index: /branches/pap_branch_051214/psModules/src/imsubtract/Makefile.am =================================================================== --- /branches/pap_branch_051214/psModules/src/imsubtract/Makefile.am (revision 5794) +++ /branches/pap_branch_051214/psModules/src/imsubtract/Makefile.am (revision 5795) @@ -4,10 +4,10 @@ libpsmoduleimsubtract_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmoduleimsubtract_la_SOURCES = pmImageSubtract.c \ - pmSubtractBias.c \ - pmSubtractSky.c + pmSubtractBias.c +# pmSubtractSky.c psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ pmImageSubtract.h \ - pmSubtractBias.h \ - pmSubtractSky.h + pmSubtractBias.h +# pmSubtractSky.h Index: /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.c =================================================================== --- /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.c (revision 5794) +++ /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.c (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmSubtractBias.c * @@ -6,28 +11,22 @@ * @author GLG, MHPCC * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-19 00:55:18 $ + * @version $Revision: 1.6.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * */ -/*****************************************************************************/ -/* INCLUDE FILES */ -/*****************************************************************************/ -#include -#include -#include -#include "pslib.h" + #if HAVE_CONFIG_H #include #endif + #include "pmSubtractBias.h" -/*****************************************************************************/ -/* DEFINE STATEMENTS */ -/*****************************************************************************/ +#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 +#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 + // XXX: put these in psConstants.h -void PS_POLY1D_PRINT( - psPolynomial1D *poly) +void PS_POLY1D_PRINT(psPolynomial1D *poly) { printf("-------------- PS_POLY1D_PRINT() --------------\n"); @@ -57,42 +56,35 @@ }\ -/*****************************************************************************/ -/* TYPE DEFINITIONS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ -psS32 currentId = 0; // XXX: remove -psS32 memLeaks = 0; // XXX: remove -//PRINT_MEMLEAKS(8); XXX -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ - /****************************************************************************** -psSubtractFrame(): this routine will take as input the pmReadout for the input -image and a pmReadout for the bias image. The bias image is subtracted in -place from the input image. We assume that sizes and types are checked -elsewhere. - -XXX: Verify that the image and readout offsets are being used the right way. - -XXX: Ensure that it does the correct thing with image size. +psSubtractFrame(): this routine will take as input a readout for the input +image and a readout for the bias image. The bias image is subtracted in +place from the input image. *****************************************************************************/ -static pmReadout *SubtractFrame( - pmReadout *in, - const pmReadout *bias) +static pmReadout *SubtractFrame(pmReadout *in, + const pmReadout *bias) { - // XXX: When did the ->row0 and ->col0 offsets get coded? - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { + psS32 i; + psS32 j; + + if (bias == NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL. Returning original image.\n"); + return(in); + } + + + if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { + psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); + return(in); + } + if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { + psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); + return(in); + } + + for (i=0;iimage->numRows;i++) { + for (j=0;jimage->numCols;j++) { in->image->data.F32[i][j]-= bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0]; - if ((in->mask != NULL) && (bias->mask != NULL)) { (in->mask->data.U8[i][j])|= @@ -105,57 +97,13 @@ } - -/****************************************************************************** -psSubtractDarkFrame(): this routine will take as input the pmReadout for the -input image and a pmReadout for the dark image. The dark image is scaled and -subtracted in place from the input image. - -XXX: Verify that the image and readout offsets are being used the right way. - -XXX: Ensure that it does the correct thing with image size. -*****************************************************************************/ -static pmReadout *SubtractDarkFrame( - pmReadout *in, - const pmReadout *dark, - psF32 scale) -{ - // XXX: When did the ->row0 and ->col0 offsets get coded? - if (fabs(scale) > FLT_EPSILON) { - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]); - - if ((in->mask != NULL) && (dark->mask != NULL)) { - (in->mask->data.U8[i][j])|= - dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; - } - } - } - } else { - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]; - - if ((in->mask != NULL) && (dark->mask != NULL)) { - (in->mask->data.U8[i][j])|= - dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; - } - } - } - } - - return(in); -} - /****************************************************************************** ImageSubtractScalar(): subtract a scalar from the input image. -XXX: Is there a psLib function for this? +XXX: Use a psLib function for this. + +XXX: This should *****************************************************************************/ -static psImage *ImageSubtractScalar( - psImage *image, - psF32 scalar) +static psImage *ImageSubtractScalar(psImage *image, + psF32 scalar) { for (psS32 i=0;inumRows;i++) { @@ -221,5 +169,5 @@ if (numOptions == 0) { - psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n"); + psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); } if (numOptions != 1) { @@ -230,95 +178,5 @@ } -/****************************************************************************** -Polynomial1DCopy(): This private function copies the members of the existing -psPolynomial1D "in" into the existing psPolynomial1D "out". The previous -members of the existing psPolynomial1D "out" are psFree'ed. - *****************************************************************************/ -static psBool Polynomial1DCopy( - psPolynomial1D *out, - psPolynomial1D *in) -{ - psFree(out->coeff); - psFree(out->coeffErr); - psFree(out->mask); - - out->type = in->type; - out->nX = in->nX; - - out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->coeff[i] = in->coeff[i]; - } - - out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->coeffErr[i] = in->coeffErr[i]; - } - - out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->mask[i] = in->mask[i]; - } - - return(true); -} - -/****************************************************************************** -Polynomial1DDup(): This private function duplicates and then returns the input -psPolynomial1D "in". - *****************************************************************************/ -static psPolynomial1D *Polynomial1DDup( - psPolynomial1D *in) -{ - psPolynomial1D *out = psPolynomial1DAlloc(in->nX, in->type); - Polynomial1DCopy(out, in); - return(out); -} - - -/****************************************************************************** -SplineCopy(): This private function copies the members of the existing -psSpline in into the existing psSpline out. - *****************************************************************************/ -static psBool SplineCopy( - psSpline1D *out, - psSpline1D *in) -{ - PS_ASSERT_PTR_NON_NULL(out, false); - PS_ASSERT_PTR_NON_NULL(in, false); - - for (psS32 i = 0 ; i < out->n ; i++) { - psFree(out->spline[i]); - } - psFree(out->spline); - psFree(out->knots); - psFree(out->p_psDeriv2); - - out->n = in->n; - out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *)); - for (psS32 i = 0 ; i < in->n ; i++) { - out->spline[i] = Polynomial1DDup(in->spline[i]); - } - - // XXX: use psVectorCopy if they get it working. - out->knots = psVectorAlloc(in->knots->n, in->knots->type.type); - for (psS32 i = 0 ; i < in->knots->n ; i++) { - out->knots->data.F32[i] = in->knots->data.F32[i]; - } - /* - out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type); - */ - - out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->n + 1) ; i++) { - out->p_psDeriv2[i] = in->p_psDeriv2[i]; - } - - return(true); -} + /****************************************************************************** @@ -328,14 +186,16 @@ PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. PM_FIT_SPLINE: fit splines to the input vector data. -The resulting spline or polynomial is set in the fitSpec argument. +XXX: Doesn't it make more sense to do polynomial interpolation on a few +elements of the input vector, rather than fit a polynomial to the entire +vector? *****************************************************************************/ -static psVector *ScaleOverscanVector( - psVector *overscanVector, - psS32 n, - void *fitSpec, - pmFit fit) +static psVector *ScaleOverscanVector(psVector *overscanVector, + psS32 n, + void *fitSpec, + pmFit fit) { psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); + // PS_VECTOR_PRINT_F32(overscanVector); if (NULL == overscanVector) { @@ -350,16 +210,20 @@ // if (n == overscanVector->n) { - return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32)); - } + for (psS32 i = 0 ; i < n ; i++) { + newVec->data.F32[i] = overscanVector->data.F32[i]; + } + return(newVec); + } + psPolynomial1D *myPoly; + psSpline1D *mySpline; psF32 x; - + psS32 i; if (fit == PM_FIT_POLYNOMIAL) { // Fit a polynomial to the old overscan vector. - psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec; + myPoly = (psPolynomial1D *) fitSpec; PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY1D(myPoly, NULL); myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n"); + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); return(NULL); } @@ -368,9 +232,17 @@ // of the old vector, use the fitted polynomial to determine the // interpolated value at that point, and set the new vector. - for (psS32 i=0;in) / ((psF32) n); newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); } } else if (fit == PM_FIT_SPLINE) { + psS32 mustFreeSpline = 0; + // Fit a spline to the old overscan vector. + mySpline = (psSpline1D *) fitSpec; + // XXX: Does it make any sense to have a psSpline argument? + if (mySpline == NULL) { + mustFreeSpline = 1; + } + // // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, @@ -378,26 +250,26 @@ // properly when doing the spline eval. // - psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector); + // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + mySpline = psVectorFitSpline1D(NULL, overscanVector); if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n"); + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); return(NULL); } + // PS_PRINT_SPLINE(mySpline); // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted spline to determine the + // of the old vector, use the fitted polynomial to determine the // interpolated value at that point, and set the new vector. - for (psS32 i=0;in - 1] x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); newVec->data.F32[i] = psSpline1DEval(mySpline, x); } - - psSpline1D *ptrSpline = (psSpline1D *) fitSpec; - if (ptrSpline != NULL) { - // Copy the resulting spline fit into ptrSpline. - PS_ASSERT_SPLINE(ptrSpline, NULL); - SplineCopy(ptrSpline, mySpline); - } - psFree(mySpline); + if (mustFreeSpline ==1) { + psFree(mySpline); + } + // PS_VECTOR_PRINT_F32(newVec); + + } else { psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); @@ -412,880 +284,405 @@ /****************************************************************************** +XXX: The SDRS does not specify type support. F32 is implemented here. *****************************************************************************/ -static psS32 GetOverscanSize( - psImage *inImg, - pmOverscanAxis overScanAxis) -{ - if (overScanAxis == PM_OVERSCAN_ROWS) { - return(inImg->numCols); - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - return(inImg->numRows); - } else if (overScanAxis == PM_OVERSCAN_ALL) { - return(1); - } - return(0); -} - -/****************************************************************************** -GetOverscanAxis(in) this private routine determines the appropiate overscan -axis from the parent cell metadata. - -XXX: Verify the READDIR corresponds with my overscan axis. - *****************************************************************************/ -static pmOverscanAxis GetOverscanAxis(pmReadout *in) -{ - psBool rc; - if ((in->parent != NULL) && (in->parent->concepts)) { - psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR"); - if (rc == true) { - if (dir == 1) { - return(PM_OVERSCAN_ROWS); - } else if (dir == 2) { - return(PM_OVERSCAN_COLUMNS); - } else if (dir == 3) { - return(PM_OVERSCAN_ALL); - } - } - } - - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata. Setting overscan axis to PM_OVERSCAN_NONE.\n"); - return(PM_OVERSCAN_NONE); -} - -/****************************************************************************** -psListLength(list): determine the length of a psList. - -XXX: Put this elsewhere. - *****************************************************************************/ -static psS32 psListLength( - psList *list) -{ - psS32 length = 0; - psListElem *tmpElem = (psListElem *) list->head; - while (NULL != tmpElem) { - tmpElem = tmpElem->next; - length++; - } - return(length); -} - -/****************************************************************************** -Note: this isn't needed anymore as of psModule SDRS 12-09. - *****************************************************************************/ -static psBool OverscanReducePixel( - psImage *in, - psList *bias, - psStats *myStats) -{ - PS_ASSERT_PTR_NON_NULL(in, NULL); - PS_ASSERT_PTR_NON_NULL(bias, NULL); - PS_ASSERT_PTR_NON_NULL(bias->head, NULL); - PS_ASSERT_PTR_NON_NULL(myStats, NULL); - - // Allocate a psVector with one element per overscan image. - psS32 numOverscanImages = psListLength(bias); - psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32); - psListElem *tmpOverscan = (psListElem *) bias->head; - psS32 i = 0; - psF64 statValue; - // - // We loop through each overscan image, calculating the specified - // statistic on that image. - // - while (NULL != tmpOverscan) { - psImage *myOverscanImage = (psImage *) tmpOverscan->data; - - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); - myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); - if (myStats == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - statsAll->data.F32[i] = statValue; - i++; - tmpOverscan = tmpOverscan->next; - } - - // - // We reduce the individual stats for each overscan image to - // a single psF32. - // - myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0); - if (myStats == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - - // - // Subtract the result and return. - // - ImageSubtractScalar(in, statValue); - psFree(statsAll); - return(in); -} - -/****************************************************************************** -ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces -a single psImage to a column by combining all pixels from each row into a -single pixel via requested statistic in myStats. - *****************************************************************************/ -static psVector *ReduceOverscanImageToCol( - psImage *overscanImage, - psStats *myStats) -{ - psF64 statValue; - psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); - psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); - - // - // For each row, we store all pixels in that row into a temporary psVector, - // then we run psVectorStats() on that vector. - // - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - tmpRow->data.F32[j] = overscanImage->data.F32[i][j]; - } - - psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); - return(NULL); - } - - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); - return(NULL); - } - - tmpCol->data.F32[i] = (psF32) statValue; - } - psFree(tmpRow); - - return(tmpCol); -} - -/****************************************************************************** -ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces -a single psImage to a row by combining all pixels from each column into a -single pixel via requested statistic in myStats. - *****************************************************************************/ -static psVector *ReduceOverscanImageToRow( - psImage *overscanImage, - psStats *myStats) -{ - psF64 statValue; - psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); - psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); - - // - // For each column, we store all pixels in that column into a temporary psVector, - // then we run psVectorStats() on that vector. - // - for (psS32 i=0;inumCols;i++) { - for (psS32 j=0;jnumRows;j++) { - tmpCol->data.F32[j] = overscanImage->data.F32[j][i]; - } - - psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); - return(NULL); - } - - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); - return(NULL); - } - - tmpRow->data.F32[i] = (psF32) statValue; - } - psFree(tmpCol); - - return(tmpRow); -} - -/****************************************************************************** -OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of -overscan images (in bias) and reduces them to a single psVector via the -specified psStats struct. The vector is then scaled to the length or the -row/column in inImg. - *****************************************************************************/ -static psVector* OverscanReduce( - psImage *inImg, - pmOverscanAxis overScanAxis, - psList *bias, - void *fitSpec, - pmFit fit, - psStats *myStats) -{ - if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) { - psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n"); - return(NULL); - } - PS_ASSERT_PTR_NON_NULL(inImg, NULL); - PS_ASSERT_PTR_NON_NULL(bias, NULL); - PS_ASSERT_PTR_NON_NULL(bias->head, NULL); - PS_ASSERT_PTR_NON_NULL(myStats, NULL); - // - // Allocate a psVector for the output of this routine. - // - psS32 vecSize = GetOverscanSize(inImg, overScanAxis); - psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32); - - // - // Allocate an array of psVectors with one psVector per element of the - // final oversan column vector. These psVectors will be used with - // psStats to reduce the multiple elements from each overscan column - // vector to a single final column vector. - // - psS32 numOverscanImages = psListLength(bias); - psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *)); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - overscanVectors[i] = NULL; - } - - // - // We iterate through the list of overscan images. For each image, - // we reduce it to a single column or row. Save the overscan vector - // in overscanVectors[]. - // - psListElem *tmpOverscan = (psListElem *) bias->head; - psS32 overscanID = 0; - while (tmpOverscan != NULL) { - psImage *tmpOverscanImage = (psImage *) tmpOverscan->data; - if (overScanAxis == PM_OVERSCAN_ROWS) { - overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats); - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats); - } - - tmpOverscan = tmpOverscan->next; - overscanID++; - } - - // - // For each overscan vector, if necessary, we scale that column or - // row to vecSize. Note: we should have already ensured that the - // fit is poly or spline. - // - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psVector *tmpOverscanVector = overscanVectors[i]; - - if (tmpOverscanVector->n != vecSize) { - overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit); - if (overscanVectors[i] == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpOverscanVector); - return(NULL); - } - psFree(tmpOverscanVector); - } - } - - // - // We collect all elements in the overscan vectors for the various - // overscan images into a single psVector (tmpVec). Then we call - // psStats on that vector to determine the final values for the - // overscan vector. - // - psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32); - psF64 statValue; - for (psS32 i = 0 ; i < vecSize ; i++) { - // Collect the i-th elements from each overscan vector into a single vector. - for (psS32 j = 0 ; j < numOverscanImages ; j++) { - tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i]; - } - - if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpVec); - return(NULL); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpVec); - return(NULL); - } - - overscanVector->data.F32[i] = (psF32) statValue; - } - - // - // We're done. Free the intermediate overscan vectors. - // - psFree(tmpVec); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - - // - // Return the computed overscanVector - // - return(overscanVector); -} - -/****************************************************************************** -RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine -takes groups of nBinOrig elements in the input vector, combines them into a -single pixel via myStats and psVectorStats(), and then outputs a vector of -those pixels. - *****************************************************************************/ -static psS32 RebinOverscanVector( - psVector *overscanVector, - psS32 nBinOrig, - psStats *myStats) -{ - psF64 statValue; - psS32 nBin; - if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { - psS32 numBins = 1+((overscanVector->n)/nBinOrig); - psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32); - psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); - - for (psS32 i=0;in > ((i*nBinOrig)+j)) { - binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; - } else { - // XXX: we get here if nBinOrig does not evenly divide - // the overscanVector vector. This is the last bin. Should - // we change the binVec->n to acknowledge that? - binVec->n = j; - } - } - psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(-1); - } - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(-1); - } - myBin->data.F32[i] = statValue; - } - - // Change the effective size of overscanVector. - overscanVector->n = numBins; - for (psS32 i=0;idata.F32[i] = myBin->data.F32[i]; - } - psFree(binVec); - psFree(myBin); - nBin = nBinOrig; - } else { - nBin = 1; - } - - return(nBin); -} - -/****************************************************************************** -FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit, -nBin): this private routine fits a psPolynomial or psSpline to the overscan -vector. It then creates a new vector, with a size determined by the input -image, evaluates the psPolynomial or psSpline at each element in that vector, -then returns that vector. - *****************************************************************************/ -static psVector *FitOverscanVectorAndUnbin( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis, - void *fitSpec, - pmFit fit, - psS32 nBin) -{ - psPolynomial1D* myPoly = NULL; - psSpline1D *mySpline = NULL; - // - // Fit a polynomial or spline to the overscan vector. - // - if (fit == PM_FIT_POLYNOMIAL) { - myPoly = (psPolynomial1D *) fitSpec; - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY1D(myPoly, NULL); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); - if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector. Returning NULL.\n"); - return(NULL); - } - } else if (fit == PM_FIT_SPLINE) { - mySpline = psVectorFitSpline1D(NULL, overscanVector); - if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning NULL.\n"); - return(NULL); - } - if (fitSpec != NULL) { - // Copy the resulting spline fit into fitSpec. - psSpline1D *ptrSpline = (psSpline1D *) fitSpec; - PS_ASSERT_SPLINE(ptrSpline, NULL); - SplineCopy(ptrSpline, mySpline); - } - } - - // - // Evaluate the poly/spline at each pixel in the overscan row/column. - // - psS32 vecSize = GetOverscanSize(inImg, overScanAxis); - psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); - if ((nBin > 1) && (nBin < overscanVector->n)) { - for (psS32 i = 0 ; i < vecSize ; i++) { - if (fit == PM_FIT_POLYNOMIAL) { - newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); - } else if (fit == PM_FIT_SPLINE) { - newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); - } - } - } else { - for (psS32 i = 0 ; i < vecSize ; i++) { - if (fit == PM_FIT_POLYNOMIAL) { - newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i); - } else if (fit == PM_FIT_SPLINE) { - newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i); - } - } - } - - psFree(mySpline); - psFree(overscanVector); - return(newVec); -} - - - -/****************************************************************************** -UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin): this private -routine takes a psVector overscanVector that was previously binned by a factor -of nBin, and then expands it to its original size, duplicated elements nBin -times for each element in the input vector overscanVector. - *****************************************************************************/ -static psVector *UnbinOverscanVector( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis, - psS32 nBin) -{ - psS32 vecSize; - - if (overScanAxis == PM_OVERSCAN_ROWS) { - vecSize = inImg->numCols; - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - vecSize = inImg->numRows; - } - - psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); - for (psS32 i = 0 ; i < vecSize ; i++) { - newVec->data.F32[i] = overscanVector->data.F32[i/nBin]; - } - - psFree(overscanVector); - return(newVec); -} - - -/****************************************************************************** -SubtractVectorFromImage(inImg, overscanVector, overScanAxis): this private -routine subtracts the overscanVector column-wise or row-wise from inImg. - *****************************************************************************/ -static psImage *SubtractVectorFromImage( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis) -{ - // - // Subtract overscan vector row-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_ROWS) { - for (psS32 i=0;inumCols;i++) { - for (psS32 j=0;jnumRows;j++) { - inImg->data.F32[j][i]-= overscanVector->data.F32[i]; - } - } - } - - // - // Subtract overscan vector column-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - inImg->data.F32[i][j]-= overscanVector->data.F32[i]; - } - } - } - - return(inImg); -} - - - -typedef enum { - PM_ERROR_NO_SUBTRACTION, - PM_WARNING_NO_SUBTRACTION, - PM_ERROR_NO_BIAS_SUBTRACT, - PM_WARNING_NO_BIAS_SUBTRACT, - PM_ERROR_NO_DARK_SUBTRACT, - PM_WARNING_NO_DARK_SUBTRACT, - PM_OKAY -} pmSubtractBiasAssertStatus; -/****************************************************************************** -AssertCodeOverscan(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for overscan subtraction. - *****************************************************************************/ -pmSubtractBiasAssertStatus AssertCodeOverscan( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - - PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - - if (overscan == true) { - pmOverscanAxis overScanAxis = GetOverscanAxis(in); - PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION); - // - // Check the type, size of each bias image. - // - psListElem *tmpOverscan = (psListElem *) in->bias->head; - psS32 numOverscans = 0; - while (NULL != tmpOverscan) { - numOverscans++; - psImage *myOverscanImage = (psImage *) tmpOverscan->data; - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); - // XXX: Get this right with the rows and columns. - if (overScanAxis == PM_OVERSCAN_ROWS) { - if (myOverscanImage->numRows != in->image->numRows) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n", - numOverscans, myOverscanImage->numCols, in->image->numRows); - if (fit == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - if (myOverscanImage->numCols != in->image->numCols) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n", - numOverscans, myOverscanImage->numCols, in->image->numCols); - if (fit == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - } else if (overScanAxis != PM_OVERSCAN_ALL) { - psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - tmpOverscan = tmpOverscan->next; - } - } else { - if (fit != PM_FIT_NONE) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n"); - return(PM_WARNING_NO_SUBTRACTION); - } - } - - // XXX: I do not like the following spec since it's useless to specify - // a psSpline as the fitSpec. - if (0) { - if ((fitSpec == NULL) && - ((fit != PM_FIT_NONE) || (overscan == true))) { - psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - - return(PM_OKAY); -} - -/****************************************************************************** -AssertCodeBias(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for bias subtraction. - *****************************************************************************/ -static pmSubtractBiasAssertStatus AssertCodeBias( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); - return(PM_ERROR_NO_BIAS_SUBTRACT); - } - if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); - return(PM_ERROR_NO_BIAS_SUBTRACT); - } - - if (bias != NULL) { - PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT); - PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); - } - return(PM_OKAY); -} - -/****************************************************************************** -AssertCodeDark(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for dark subtraction. - *****************************************************************************/ -pmSubtractBiasAssertStatus AssertCodeDark( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) { - psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows. Returning in image\n"); - return(PM_ERROR_NO_DARK_SUBTRACT); - } - if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) { - psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns. Returning in image\n"); - return(PM_ERROR_NO_DARK_SUBTRACT); - } - - if (dark != NULL) { - PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT); - PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); - } - return(PM_OKAY); -} - -/****************************************************************************** -p_psDetermineTrimmedImage(): global routine: determines the region of the -input pmReadout which will be operated on by the various detrend modules. It -does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout. - -Use it this way: - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - -XXX: Create a pmUtils.c file and put this routine there. - *****************************************************************************/ -psImage *p_psDetermineTrimmedImage(pmReadout *in) -{ - if ((in->parent == NULL) || (in->parent->concepts == NULL)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n"); - return(in->image); - } - - psBool rc = false; - psImage *trimmedImg = NULL; - psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts, - "CELL.TRIMSEC"); - if (rc == false) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n"); - trimmedImg = in->image; - } else { - trimmedImg = psImageSubset(in->image, *trimRegion); - } - - return(trimmedImg); -} - - - - - - - - - - - - - - - - - - -/****************************************************************************** -pmSubtractBias(....): see SDRS for complete specification. - -XXX: Code and assert type support: U16, S32, F32. -XXX: Add trace messages. - *****************************************************************************/ -pmReadout *pmSubtractBias( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBin, - const pmReadout *bias, - const pmReadout *dark) +pmReadout *pmSubtractBias(pmReadout *in, + void *fitSpec, + const psList *overscans, + pmOverscanAxis overScanAxis, + psStats *stat, + psS32 nBinOrig, + pmFit fit, + const pmReadout *bias) { psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, "---- pmSubtractBias() begin ----\n"); - // - // Check input parameters, generate warnings and errors. - // - if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - return(in); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - - // - // Subtract overscan frames if necessary. - // - if (overscan == true) { - pmOverscanAxis overScanAxis = GetOverscanAxis(in); + PS_ASSERT_READOUT_NON_NULL(in, NULL); + PS_ASSERT_READOUT_NON_EMPTY(in, NULL); + PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); + + // + // If the overscans != NULL, then check the type of each image. + // + if (overscans != NULL) { + psListElem *tmpOverscan = (psListElem *) overscans->head; + while (NULL != tmpOverscan) { + psImage *myOverscanImage = (psImage *) tmpOverscan->data; + PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); + tmpOverscan = tmpOverscan->next; + } + } + + if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) { + psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE). Returning in image\n"); + return(in); + } + + // Check for an unallowable pmFit. + if ((fit != PM_OVERSCAN_NONE) && + (fit != PM_OVERSCAN_ROWS) && + (fit != PM_OVERSCAN_COLUMNS) && + (fit != PM_OVERSCAN_ALL)) { + psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d). Returning in image.\n", fit); + return(in); + } + // Check for an unallowable pmOverscanAxis. + if ((overScanAxis != PM_OVERSCAN_NONE) && + (overScanAxis != PM_OVERSCAN_ROWS) && + (overScanAxis != PM_OVERSCAN_COLUMNS) && + (overScanAxis != PM_OVERSCAN_ALL)) { + psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d). Returning in image.\n", overScanAxis); + return(in); + } + psS32 i; + psS32 j; + psS32 numBins = 0; + static psVector *overscanVector = NULL; + psVector *tmpRow = NULL; + psVector *tmpCol = NULL; + psVector *myBin = NULL; + psVector *binVec = NULL; + psListElem *tmpOverscan = NULL; + double statValue; + psImage *myOverscanImage = NULL; + psPolynomial1D *myPoly = NULL; + psSpline1D *mySpline = NULL; + psS32 nBin; + + // + // Create a static stats data structure and determine the highest + // priority stats option. + // + static psStats *myStats = NULL; + if (myStats == NULL) { + myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + p_psMemSetPersistent(myStats, true); + } + if (stat != NULL) { + myStats->options = GenNewStatOptions(stat); + } + + + if (overScanAxis == PM_OVERSCAN_NONE) { + if (fit != PM_FIT_NONE) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE. Proceeding to full fram subtraction.\n"); + } + + if (overscans != NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL. Proceeding to full fram subtraction.\n"); + } + return(SubtractFrame(in, bias)); + } + + if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE. Proceeding with the rest of the module.\n"); + } + + + // + // We subtract each overscan region from the image data. + // If we get here we know that overscans != NULL. + // + + if (overScanAxis == PM_OVERSCAN_ALL) { + tmpOverscan = (psListElem *) overscans->head; + while (NULL != tmpOverscan) { + myOverscanImage = (psImage *) tmpOverscan->data; + + PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); + psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(myStats, &statValue)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); + return(in); + } + ImageSubtractScalar(in->image, statValue); + + tmpOverscan = tmpOverscan->next; + } + return(in); + } + + // This check is redundant with above code. + if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) { + psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis); + return(in); + } + + tmpOverscan = (psListElem *) overscans->head; + while (NULL != tmpOverscan) { + // PS_IMAGE_PRINT_F32_HIDEF(in->image); + myOverscanImage = (psImage *) tmpOverscan->data; + + if (overScanAxis == PM_OVERSCAN_ROWS) { + if (myOverscanImage->numCols != (in->image)->numCols) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n", + myOverscanImage->numCols, in->image->numCols); + } + + // We create a row vector and subtract this vector from image. + // XXX: Is there a better way to extract a psVector from a psImage without + // having to copy every element in that vector? + overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); + for (i=0;in;i++) { + overscanVector->data.F32[i] = 0.0; + } + tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); + + // For each column of the input image, loop through every row, + // collect the pixel in that row, then performed the specified + // statistical op on those pixels. Store this in overscanVector. + for (i=0;inumCols;i++) { + for (j=0;jnumRows;j++) { + tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i]; + } + psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(rc, &statValue)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); + return(in); + } + overscanVector->data.F32[i] = statValue; + } + psFree(tmpRow); + + // Scale the overscan vector to the size of the input image. + if (overscanVector->n != in->image->numCols) { + if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { + psVector *newVec = ScaleOverscanVector(overscanVector, + in->image->numCols, + fitSpec, fit); + if (newVec == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); + return(in); + } + psFree(overscanVector); + overscanVector = newVec; + } else { + psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); + psFree(overscanVector); + return(in); + } + } + } + + if (overScanAxis == PM_OVERSCAN_COLUMNS) { + if (myOverscanImage->numRows != (in->image)->numRows) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n", + myOverscanImage->numRows, in->image->numRows); + } + + // We create a column vector and subtract this vector from image. + overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); + for (i=0;in;i++) { + overscanVector->data.F32[i] = 0.0; + } + tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); + + // For each row of the input image, loop through every column, + // collect the pixel in that row, then performed the specified + // statistical op on those pixels. Store this in overscanVector. + for (i=0;inumRows;i++) { + for (j=0;jnumCols;j++) { + tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j]; + } + psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(rc, &statValue)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); + return(in); + } + overscanVector->data.F32[i] = statValue; + } + psFree(tmpCol); + + // Scale the overscan vector to the size of the input image. + if (overscanVector->n != in->image->numRows) { + if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { + psVector *newVec = ScaleOverscanVector(overscanVector, + in->image->numRows, + fitSpec, fit); + if (newVec == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); + return(in); + } + psFree(overscanVector); + overscanVector = newVec; + } else { + psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); + psFree(overscanVector); + return(in); + } + } + } + // - // Create a psStats data structure and determine the highest - // priority stats option. + // Re-bin the overscan vector (change its length). // - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - if (stat != NULL) { - myStats->options = GenNewStatOptions(stat); - } - + // Only if nBinOrig > 1. + if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { + numBins = 1+((overscanVector->n)/nBinOrig); + myBin = psVectorAlloc(numBins, PS_TYPE_F32); + binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); + + for (i=0;in > ((i*nBinOrig)+j)) { + binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; + } else { + // XXX: we get here if nBinOrig does not evenly divide + // the overscanVector vector. This is the last bin. Should + // we change the binVec->n to acknowledge that? + binVec->n = j; + } + } + psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); + if (rc == NULL) { + psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); + return(in); + } + if (false == p_psGetStatValue(rc, &statValue)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); + return(in); + } + myBin->data.F32[i] = statValue; + } + + // Change the effective size of overscanVector. + overscanVector->n = numBins; + for (i=0;idata.F32[i] = myBin->data.F32[i]; + } + psFree(binVec); + psFree(myBin); + nBin = nBinOrig; + } else { + nBin = 1; + } + + // At this point the number of data points in overscanVector should be + // equal to the number of rows/columns (whatever is appropriate) in the + // image divided by numBins. // - // Reduce overscan images to a single pixel, then subtract. - // This code is no longer required as of SDRS 12-09. + + // - if (overScanAxis == PM_OVERSCAN_ALL) { - if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) { - return(in); - } - psFree(myStats); + // This doesn't seem right. The only way to do a spline fit is if, + // by SDRS requirements, fitSpec is not-NULL> But in order for it + // to be non-NULL, someone must have called psSpline1DAlloc() with + // the min, max, and number of splines. + // + if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) { + // + // Fit a polynomial or spline to the overscan vector. + // + if (fit == PM_FIT_POLYNOMIAL) { + myPoly = (psPolynomial1D *) fitSpec; + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); + if (myPoly == NULL) { + psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector. Returning in image.\n"); + psFree(overscanVector); + return(in); + } + } else if (fit == PM_FIT_SPLINE) { + // XXX: This makes no sense + // XXX: must free mySpline? + mySpline = (psSpline1D *) fitSpec; + mySpline = psVectorFitSpline1D(NULL, overscanVector); + if (mySpline == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning in image.\n"); + psFree(overscanVector); + return(in); + } + } + + // + // Subtract fitted overscan vector row-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_ROWS) { + for (i=0;i<(in->image)->numCols;i++) { + psF32 tmpF32 = 0.0; + if (fit == PM_FIT_POLYNOMIAL) { + tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); + } else if (fit == PM_FIT_SPLINE) { + tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); + } + for (j=0;j<(in->image)->numRows;j++) { + (in->image)->data.F32[j][i]-= tmpF32; + } + } + } + + // + // Subtract fitted overscan vector column-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_COLUMNS) { + for (i=0;i<(in->image)->numRows;i++) { + psF32 tmpF32 = 0.0; + if (fit == PM_FIT_POLYNOMIAL) { + tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); + } else if (fit == PM_FIT_SPLINE) { + tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); + } + + for (j=0;j<(in->image)->numCols;j++) { + (in->image)->data.F32[i][j]-= tmpF32; + } + } + } } else { // - // Reduce the overscan images to a single overscan vector. - // - psVector *overscanVector = OverscanReduce(in->image, overScanAxis, - in->bias, fitSpec, - fit, myStats); - if (overscanVector == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector. Returning in image\n"); - psFree(myStats); - return(in); - } - - // - // Rebin the overscan vector if necessary. - // - psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats); - if (newBin < 0) { - psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector. Returning in image\n"); - psFree(myStats); - return(in); - } - - // - // If necessary, fit a psPolynomial or psSpline to the overscan vector. - // Then, unbin the overscan vector to appropriate length for the in image. - // - if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { - overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin); - if (overscanVector == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector. Returning in image\n"); - psFree(myStats); - return(in); - } - } else { - overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin); - } - - // - // Subtract the overscan vector from the input image. - // - SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis); - psFree(myStats); - psFree(overscanVector); - } - } - - // - // Perform bias subtraction if necessary. - // - if (bias != NULL) { - if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - SubtractFrame(in, bias); - } - } - - // - // Perform dark subtraction if necessary. - // - if (dark != NULL) { - if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - psBool rc; - psF32 scale = 0.0; - if (in->parent != NULL) { - scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME"); - if (rc == false) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n"); - } - } - SubtractDarkFrame(in, dark, scale); - } - } - - // - // All done. - // + // If we get here, then no polynomials were fit to the overscan + // vector. We simply subtract it, taking into account binning, + // from the image. + // + + // + // Subtract overscan vector row-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_ROWS) { + for (i=0;i<(in->image)->numCols;i++) { + for (j=0;j<(in->image)->numRows;j++) { + (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin]; + } + } + } + + // + // Subtract overscan vector column-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_COLUMNS) { + for (i=0;i<(in->image)->numRows;i++) { + for (j=0;j<(in->image)->numCols;j++) { + (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin]; + } + } + } + } + + psFree(overscanVector); + + tmpOverscan = tmpOverscan->next; + } + psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, "---- pmSubtractBias() exit ----\n"); + + if (bias != NULL) { + return(SubtractFrame(in, bias)); + } return(in); } Index: /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.h =================================================================== --- /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.h (revision 5794) +++ /branches/pap_branch_051214/psModules/src/imsubtract/pmSubtractBias.h (revision 5795) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmSubtractBias.h * @@ -6,6 +11,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-23 23:54:30 $ + * @version $Revision: 1.4.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-17 03:18:39 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -23,4 +28,5 @@ #include #include "pslib.h" + #include "pmAstrometry.h" @@ -38,25 +44,12 @@ } pmFit; -pmReadout *pmSubtractBias( - pmReadout *in, ///< The input pmReadout image - void *fitSpec, ///< A polynomial or spline, defining the fit type. - pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE - bool overscan, - psStats *stat, ///< The statistic to be used in combining overscan data - int nBin, ///< The amount of binning to be done image pixels. - const pmReadout *bias, ///< A possibly NULL bias pmReadout which is to be subtracted - const pmReadout *dark ///< A possibly NULL bias pmReadout which is to be subtracted -); - -/****************************************************************************** -DetermineTrimmedImageg() This private routine determines the region of the -input pmReadout which will be operated on by the various detrend modules. It -does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout. - -XXX: Consider making a pmUtils.c file and put this routine there. - *****************************************************************************/ -psImage *p_psDetermineTrimmedImage( - pmReadout *in -); +pmReadout *pmSubtractBias(pmReadout *in, ///< The input pmReadout image + void *fitSpec, ///< A polynomial or spline, defining the fit type. + const psList *overscans, ///< A psList of overscan images + pmOverscanAxis overScanAxis, ///< Defines how overscans are applied + psStats *stat, ///< The statistic to be used in combining overscan data + int nBin, ///< The amount of binning to be done image pixels. + pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE + const pmReadout *bias); ///< A possibly NULL bias pmReadout which is to be subtracted #endif