Index: /trunk/archive/scripts/src/phase2/pmFlatField.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatField.c (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmFlatField.c (revision 5107) @@ -0,0 +1,177 @@ +/** @file pmFlatField.c + * + * @brief Given an input image and a flat field image, pmFlatField shall divide the input image by the flat + * field image. + * + * The input image, in, and the flat field image, flat, need not be the same size, since the input image may + * already have been trimmed (following overscan subtraction), but the function shall use the offsets in the + * image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct pixel on the flat + * field. In the event that the flat image is too small (i.e., pixels on the input image refer to pixels + * outside the range of the flat image), the function shall generate an error. Pixels which are negative or + * zero in the flat shall be masked in the input image with the value PM_MASK_FLAT. Negative pixels in the + * flat may be set to zero so that they are treated identically to zeroes. Any pixels masked in the flat + * shall be masked with corresponding values in the output. The function shall not normalize the flat; this + * responsibility is left to the caller. This function is basically equivalent to a divide (with psImageOp), + * but with care for the region that is divided, checking for negative pixels, and copying of the mask from + * the flat to the output. + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include +#include "pslib.h" + +#include "pmFPA.h" +#include "pmFlatField.h" +#include "pmMaskBadPixels.h" +#include "pmFlatFieldErrors.h" + + +bool pmFlatField(pmReadout *in, const pmReadout *flat) +{ + int i = 0; + int j = 0; + int totOffCol = 0; + int totOffRow = 0; + psElemType inType; + psElemType flatType; + psElemType maskType; + psImage *inImage = NULL; + psImage *inMask = NULL; + psImage *flatImage = NULL; + + + // 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; + 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; + } 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 flat image + if(inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE, + inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols); + return false; + } else if(inMask->numRows>flatImage->numRows || inMask->numCols > flatImage->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE, + inMask->numRows, inMask->numCols, flatImage->numRows, flatImage->numCols); + return false; + } + + // Determine total offset based on image offset with chip offset + totOffCol = inImage->col0 + in->col0; + totOffRow = inImage->row0 + in->row0; + + // Check that offsets are within image limits + if(totOffRow>=flatImage->numRows || totOffCol>=flatImage->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmFlatField_OFFSET_FLAT_IMAGE, + totOffRow, totOffCol, flatImage->numRows, flatImage->numCols); + return false; + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE, + totOffRow, totOffCol, inImage->numRows, inImage->numCols); + return false; + } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE, + totOffRow, totOffCol, inMask->numRows, inMask->numCols); + return false; + } + + // Check for incorrect types + inType = inImage->type.type; + flatType = flatImage->type.type; + maskType = inMask->type.type; + if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmFlatField_TYPE_INPUT_IMAGE, + inType); + return false; + } else if(PS_IS_PSELEMTYPE_COMPLEX(flatType)) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmFlatField_TYPE_FLAT_IMAGE, + flatType); + return false; + } else if(maskType != PS_TYPE_MASK) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE, + maskType); + return false; + } else if(inType != flatType) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmFlatField_TYPE_MISMATCH, + inType, flatType); + return false; + } + + // Macro for all PS types + #define PM_FLAT_DIVISION(TYPE) \ +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 < 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 */ \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ + flatImage->data.TYPE[j][i] = 0.0; \ + } \ + } \ + } \ + 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 */ \ + inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ + } \ + } \ + } \ + break; + + switch(inType) { + PM_FLAT_DIVISION(U8); + PM_FLAT_DIVISION(U16); + PM_FLAT_DIVISION(U32); + PM_FLAT_DIVISION(U64); + PM_FLAT_DIVISION(S8); + PM_FLAT_DIVISION(S16); + PM_FLAT_DIVISION(S32); + PM_FLAT_DIVISION(S64); + PM_FLAT_DIVISION(F32); + PM_FLAT_DIVISION(F64); + default: + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmFlatField_TYPE_UNSUPPORTED, + inType); + } + + return true; +} Index: /trunk/archive/scripts/src/phase2/pmFlatField.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatField.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmFlatField.h (revision 5107) @@ -0,0 +1,40 @@ +/** @file pmFlatField.h + * + * @brief Given an input image and a flat field image, pmFlatField shall divide the input image by the flat + * field image. + * + * The input image, in, and the flat field image, flat, need not be the same size, since the input image may + * already have been trimmed (following overscan subtraction), but the function shall use the offsets in the + * image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct pixel on the flat + * field. In the event that the flat image is too small (i.e., pixels on the input image refer to pixels + * outside the range of the flat image), the function shall generate an error. Pixels which are negative or + * zero in the flat shall be masked in the input image with the value PM_MASK_FLAT. Negative pixels in the + * flat may be set to zero so that they are treated identically to zeroes. Any pixels masked in the flat + * shall be masked with corresponding values in the output. The function shall not normalize the flat; this + * responsibility is left to the caller. This function is basically equivalent to a divide (with psImageOp), + * but with care for the region that is divided, checking for negative pixels, and copying of the mask from + * the flat to the output. + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +/** Execute flat field module. + * + * Given an input image and a flat-field image, pmFlatField shall divide the input image by the flat field + * image. + * + * @return bool: True or false for success or failure + */ + +#include "pmFPA.h" + +bool pmFlatField( + pmReadout *in, ///< Redout with input image and mask + const pmReadout *flat ///< Readout with flat image +); + Index: /trunk/archive/scripts/src/phase2/pmFlatFieldErrors.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmFlatFieldErrors.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmFlatFieldErrors.h (revision 5107) @@ -0,0 +1,47 @@ +/** @file pmFlatFieldErrors.h + * + * @brief Contains the error text for the flat field module + * + * @ingroup ErrorHandling + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:59:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PM_FLATFIELD_ERRORS_H +#define PM_FLATFIELD_ERRORS_H + +/* N.B., lines between '//~Start' and '//~End' are automatic generated from + * the template following the '//~Start'. The template is used to generate + * the other lines by, for each error text in psDataManipErrors.dat, the following + * substitutions are made: + * $1 The error text macro name (first word in the psFlatFieldErrors.h lines) + * $2 The error text (rest of the line in psFlatFieldErrors.h) + * $n The order of the source line in psFlatFieldErrors.h (comments excluded) + * + * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. + */ + +#define PS_ERRORNAME_DOMAIN "psModule.src." + +//~Start #define PS_ERRORTEXT_$1 "$2" +#define PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT "Null not allowed for flat readout." +#define PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE "Null not allowed for input image." +#define PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE "Null not allowed for flat image." +#define PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE "Input image size exceeds that of flat image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE "Input image mask size exceeds that of flat image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmFlatField_OFFSET_FLAT_IMAGE "Total offset >= flat image size: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE "Total offset >= input image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE "Total offset >= input image mask: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmFlatField_TYPE_INPUT_IMAGE "Complex types not allowed for input image. Type: %d" +#define PS_ERRORTEXT_pmFlatField_TYPE_FLAT_IMAGE "Complex types not allowed for flat image. Type: %d" +#define PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE "Mask must be PS_TYPE_MASK type. Type: %d" +#define PS_ERRORTEXT_pmFlatField_TYPE_MISMATCH "Input and flat image types differ: (%d vs %d)" +#define PS_ERRORTEXT_pmFlatField_TYPE_UNSUPPORTED "Unsupported image datatype. Type: %d" +//~End + +#endif Index: /trunk/archive/scripts/src/phase2/pmMaskBadPixels.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixels.c (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmMaskBadPixels.c (revision 5107) @@ -0,0 +1,183 @@ +/** @file pmMaskBadPixels.c + * + * @brief Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a + * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that + * match the value to mask. + * + * Given an input image, in, a bad pixel mask, a corresponding value in the bad pixel mask to mask in the + * input image, maskVal, a saturation level, and a growing radius, pmMaskBadPixels shall mask in the input + * image those pixels in the bad pixel mask that match the value to mask. Note that the input image, in, is + * modified in-place. All pixels in the mask which satisfy the maskVal shall have their corresponding pixels + * masked in the input image, in. All pixels which satisfy the growVal shall have their corresponding + * pixels, along with all pixels within the grow radius masked. Pixels which have flux greater than sat shall + * also be masked, but not grown. Note that the input image, in, and the mask need not be the same size, since + * the input image may already have been trimmed (following overscan subtraction), but the function shall use + * the offsets in the image (in->x0 and in->y0) to determine the appropriate offsets to obtain the correct + * pixel on the mask. In the event that the mask image is too small (i.e., pixels on the input image refer to + * pixels outside the range of the mask image), the function shall generate an error. + + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmMaskBadPixels.h" +#include "pmMaskBadPixelsErrors.h" + +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; + + + // 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; + } + + 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(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE, + inImage->numRows, inImage->numCols, mask->numRows, mask->numCols); + return false; + } 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->numRows, mask->numCols); + return false; + } + + // Determine total offset based on image offset with chip offset + totOffCol = inImage->col0 + in->col0; + totOffRow = inImage->row0 + in->row0; + + // Check that offsets are within image limits + if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE, + totOffRow, totOffCol, mask->numRows, mask->numCols); + return false; + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE, + totOffRow, totOffCol, inImage->numRows, inImage->numCols); + return false; + } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { + psError( PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE_MASK, + totOffRow, totOffCol, inMask->numRows, inMask->numCols); + return false; + } + + // Check for incorrect types + inType = inImage->type.type; + maskType = mask->type.type; + if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmMaskBadPixels_TYPE_INPUT_IMAGE, + inType); + return false; + } else if(maskType!=PS_TYPE_MASK) { + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE, + maskType); + return false; + } + + // 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; + + // 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 false; +} Index: /trunk/archive/scripts/src/phase2/pmMaskBadPixels.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixels.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmMaskBadPixels.h (revision 5107) @@ -0,0 +1,61 @@ +/** @file pmMaskBadPixels.h + * + * @brief Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a + * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that + * match the value to mask. + * + * Given an input image, in, a bad pixel mask, a corresponding value in the bad pixel mask to mask in the + * input image, maskVal, a saturation level, and a growing radius, pmMaskBadPixels shall mask in the input + * image those pixels in the bad pixel mask that match the value to mask. Note that the input image, in, is + * modified in-place. All pixels in the mask which satisfy the maskVal shall have their corresponding pixels + * masked in the input image, in. All pixels which satisfy the growVal shall have their corresponding + * pixels, along with all pixels within the grow radius masked. Pixels which have flux greater than sat shall + * also be masked, but not grown. Note that the input image, in, and the mask need not be the same size, + * since the input image may already have been trimmed (following overscan subtraction), but the function + * shall use the offsets in the image (in->x0 and in->y0) to determine the appropriate offsets to obtain the + * correct pixel on the mask. In the event that the mask image is too small (i.e., pixels on the input image + * refer to pixels outside the range of the mask image), the function shall generate an error. + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#include "pslib.h" +#include "pmFPA.h" + +/** Mask values */ +typedef enum { + PM_MASK_TRAP = 0x0001, ///< The pixel is a charge trap. + PM_MASK_BADCOL = 0x0002, ///< The pixel is a bad column. + PM_MASK_SAT = 0x0004, ///< The pixel is saturated. + PM_MASK_FLAT = 0x0008 ///< The pixel is non-positive in the flat-field. +} pmMaskValue; + +/** Macro to find maximum of two numbers */ +#define MAX(A,B)((A)>=(B)?(A):(B)) + +/** Macro to find minimum of two numbers */ +#define MIN(A,B)((A)<=(B)?(A):(B)) + + +/** Execute bad pixels module. + * + * Given an input image, a bad pixel mask, a corresponding value in the bad pixel mask to mask, a + * saturation level, and a growing radius, mask in the input image those pixels in the bad pixel mask that + * match the value to mask. + * + * @return bool: True or false for success or failure + */ +bool pmMaskBadPixels( + pmReadout *in, ///< Readout containing input image 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. + unsigned int growVal, ///< Mask data to determine if a circurlar area should be masked. + int grow ///< Radius of mask to apply around pixel. +); + Index: /trunk/archive/scripts/src/phase2/pmMaskBadPixelsErrors.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmMaskBadPixelsErrors.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmMaskBadPixelsErrors.h (revision 5107) @@ -0,0 +1,45 @@ +/** @file pmMaskBadPixelsErrors.h + * + * @brief Contains the error text for the mask bad pixels module + * + * @ingroup ErrorHandling + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PM_FLATFIELD_ERRORS_H +#define PM_FLATFIELD_ERRORS_H + +/* N.B., lines between '//~Start' and '//~End' are automatic generated from + * the template following the '//~Start'. The template is used to generate + * the other lines by, for each error text in pmMaskBadPixelsErrors.dat, the following + * substitutions are made: + * $1 The error text macro name (first word in the pmMaskBadPixelsErrors.h lines) + * $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. + */ + +#define PS_ERRORNAME_DOMAIN "psModule.src." + +//~Start #define PS_ERRORTEXT_$1 "$2" +#define PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE "Null not allowed for mask image." +#define PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE "Null not allowed for input image." +#define PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE "Input image size exceeds that of mask image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE "Input image mask size exceeds that of mask image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE "Total offset >= mask image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE "Total offset >= input image: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE_MASK "Total offset >= input image mask: (%d, %d) vs (%d, %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_INPUT_IMAGE "Complex types not allowed for input image. Type: %d" +#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE "Mask must be PS_TYPE_MASK type. Type: %d" +#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_MISMATCH "Input and flat image types differ: (%d vs %d)" +#define PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED "Unsupported image datatype. Type: %d" +//~End + +#endif Index: /trunk/archive/scripts/src/phase2/pmNonLinear.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmNonLinear.c (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmNonLinear.c (revision 5107) @@ -0,0 +1,123 @@ +/** @file pmNonLinear.c + * + * Provides polynomial or table lookup non-linearity corrections to readouts. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: The SDR is silent about image types. Only F32 was implemented. + * + */ + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmNonLinear.h" + +/****************************************************************************** +pmNonLinearityLookup(): This routine will take an psReadout 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 +that value. + *****************************************************************************/ + +pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, + const psPolynomial1D *input1DPoly) +{ + 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_ASSERT_PTR_NON_NULL(input1DPoly, NULL); + + psS32 i; + psS32 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]); + } + } + return(inputReadout); +} + + +/****************************************************************************** +pmNonLinearityLookup(): This routine will take an psReadout image as input +and two input vectors, which constitute a lookup table. For each pixel in +the input image, that pixels value will be determined in the input vector +inFluxe, and the corresponding value in outFlux. The image pixel will then +be set to the value from outFlux. + *****************************************************************************/ +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_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); + } + 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); + + psS32 i; + psS32 j; + psS32 binNum; + psScalar x; + psS32 numPixels = 0; + 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); + } + return(inputReadout); +} Index: /trunk/archive/scripts/src/phase2/pmNonLinear.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmNonLinear.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmNonLinear.h (revision 5107) @@ -0,0 +1,27 @@ +/** @file pmNonLinear.h + * + * Provides polynomial or table lookup non-linearity corrections to readouts. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_NON_LINEAR_H) +#define PM_NON_LINEAR_H + +#include "pslib.h" +#include "pmFPA.h" + +pmReadout *pmNonLinearityPolynomial(pmReadout *in, + const psPolynomial1D *coeff); + +pmReadout *pmNonLinearityLookup(pmReadout *in, + const psVector *inFlux, + const psVector *outFlux); + +#endif Index: /trunk/archive/scripts/src/phase2/pmSubtractBias.c =================================================================== --- /trunk/archive/scripts/src/phase2/pmSubtractBias.c (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmSubtractBias.c (revision 5107) @@ -0,0 +1,646 @@ +/** @file pmSubtractBias.c + * + * This file will contain a module which will subtract the detector bias + * in place from an input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if HAVE_CONFIG_H +#include +#endif + +#include "pmSubtractBias.h" + +#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 +#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 + +/****************************************************************************** +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) +{ + 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])|= + bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0]; + } + } + } + + return(in); +} + +/****************************************************************************** +ImageSubtractScalar(): subtract a scalar from the input image. + +XXX: Use a psLib function for this. + +XXX: This should + *****************************************************************************/ +static psImage *ImageSubtractScalar(psImage *image, + psF32 scalar) +{ + for (psS32 i=0;inumRows;i++) { + for (psS32 j=0;jnumCols;j++) { + image->data.F32[i][j]-= scalar; + } + } + return(image); +} + +/****************************************************************************** +GenNewStatOptions(): this routine will take as input the options member of the +stat data structure, determine if multiple options have been specified, issue +a warning message if so, and return the highest priority option (according to +the order of the if-statements in this code). The higher priority options are +listed lower in the code. + *****************************************************************************/ +static psStatsOptions GenNewStatOptions(const psStats *stat) +{ + psS32 numOptions = 0; + psStatsOptions opt = 0; + + if (stat->options & PS_STAT_ROBUST_MODE) { + if (numOptions == 0) { + opt = PS_STAT_ROBUST_MODE; + } + numOptions++; + } + + if (stat->options & PS_STAT_ROBUST_MEDIAN) { + if (numOptions == 0) { + opt = PS_STAT_ROBUST_MEDIAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_ROBUST_MEAN) { + if (numOptions == 0) { + opt = PS_STAT_ROBUST_MEAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_CLIPPED_MEAN) { + if (numOptions == 0) { + opt = PS_STAT_CLIPPED_MEAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_SAMPLE_MEDIAN) { + if (numOptions == 0) { + opt = PS_STAT_SAMPLE_MEDIAN; + } + numOptions++; + } + + if (stat->options & PS_STAT_SAMPLE_MEAN) { + numOptions++; + opt = PS_STAT_SAMPLE_MEAN; + } + + + if (numOptions == 0) { + psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); + } + if (numOptions != 1) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.c: GenNewStatOptions(): Too many statistics options have been specified\n"); + } + return(opt); +} + +/****************************************************************************** +ScaleOverscanVector(): this routine takes as input an arbitrary vector, +creates a new vector of length n, and fills the new vector with the +interpolated values of the old vector. The type of interpolation is: + PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. + PM_FIT_SPLINE: fit splines to the input vector data. +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) +{ + psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, + "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); + + if (NULL == overscanVector) { + return(overscanVector); + } + + // Allocate the new vector. + psVector *newVec = psVectorAlloc(n, PS_TYPE_F32); + + // + // If the new vector is the same size as the old, simply copy the data. + // + if (n == overscanVector->n) { + 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. + myPoly = (psPolynomial1D *) fitSpec; + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); + if (myPoly == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); + return(NULL); + } + + // For each element of the new vector, convert the x-ordinate to that + // of the old vector, use the fitted polynomial to determine the + // interpolated value at that point, and set the new vector. + for (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; + if (mySpline == NULL) { + mustFreeSpline = 1; + } + + // + // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, + // splines enpoints will be from 0.0 to overscanVector->n-1. Must scale + // properly when doing the spline eval. + // + mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + if (mySpline == NULL) { + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); + return(NULL); + } + + // For each element of the new vector, convert the x-ordinate to that + // of the old vector, use the fitted polynomial to determine the + // interpolated value at that point, and set the new vector. + for (i=0;in - 1] + x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); + newVec->data.F32[i] = psSpline1DEval(mySpline, x); + } + if (mustFreeSpline ==1) { + psFree(mySpline); + } + + } else { + psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); + psFree(newVec); + return(NULL); + } + + psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, + "---- ScaleOverscanVector() exit ----\n"); + return(newVec); +} + +/****************************************************************************** +XXX: The SDRS does not specify type support. F32 is implemented here. + *****************************************************************************/ +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"); + 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) { + 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); + } + } + } + + // + // Re-bin the overscan vector (change its length). + // + // 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. + // + + + // + // 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) { + mySpline = (psSpline1D *) fitSpec; + mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + 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 { + // + // 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: /trunk/archive/scripts/src/phase2/pmSubtractBias.h =================================================================== --- /trunk/archive/scripts/src/phase2/pmSubtractBias.h (revision 5107) +++ /trunk/archive/scripts/src/phase2/pmSubtractBias.h (revision 5107) @@ -0,0 +1,49 @@ +/** @file pmSubtractBias.h + * + * This file will contain a module which will subtract the detector bias + * in place from an input image. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-23 02:58:30 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_SUBTRACT_BIAS_H) +#define PM_SUBTRACT_BIAS_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pmFPA.h" + +typedef enum { + PM_OVERSCAN_NONE, ///< No overscan subtraction + PM_OVERSCAN_ROWS, ///< Subtract rows + PM_OVERSCAN_COLUMNS, ///< Subtract columns + PM_OVERSCAN_ALL ///< Subtract the statistic of all pixels in overscan region +} pmOverscanAxis; + +typedef enum { + PM_FIT_NONE, ///< No fit + PM_FIT_POLYNOMIAL, ///< Fit polynomial + PM_FIT_SPLINE ///< Fit cubic splines +} pmFit; + +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