Index: /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am (revision 14720) +++ /branches/eam_branch_20070830/psLib/src/imageops/Makefile.am (revision 14721) @@ -15,4 +15,6 @@ psImageMaskOps.c \ psImageBinning.c \ + psImageMap.c \ + psImagePixelInterpolate.c \ psImageUnbin.c @@ -30,4 +32,6 @@ psImageMaskOps.h \ psImageBinning.h \ + psImageMap.h \ + psImagePixelInterpolate.h \ psImageUnbin.h Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14721) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14721) @@ -0,0 +1,110 @@ +/** @file psImageMap.c + * + * @brief Functions define a 2d coarse representation of a finer 2D field + * + * @ingroup Image + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-01 00:41:01 $ + * + * Copyright 2007 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include +#include "psError.h" +#include "psAbort.h" +#include "psAssert.h" +#include "psRegion.h" +#include "psImage.h" +#include "psImageBinning.h" +#include "psImageMap.h" + +static void psImageMapFree(psImageMap *map) { + + if (!map) return; + + psFree (map->map); + psFree (map->field); + psFree (map->stats); + psFree (map->binning); + + return; +} + +psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) { + + assert (binning); + assert (stats); + + psImageMap *map = (psImageMap*)psAlloc(sizeof(psImageMap)); + psMemSetDeallocator(map, (psFreeFunc)psImageMapFree); + + map->binning = psMemIncrRefCounter (binning); + map->field = psMemIncrRefCounter (field); + map->stats = psMemIncrRefCounter (stats); + map->map = psImageAlloc (binning->nXruff, binning->nYruff, PS_TYPE_F32); + + psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER); + + return map; +} + +// generate a psImageMap (or NULL) with the given number of superpixels in X and Y +psImageMap *psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac) { + + int nBad; + int nPoor; + + psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK); + + // accumulate the values for each map pixel + + // we can do this by accumulating a sub-vector for each cell + psArray *vectors = psArrayAlloc (map->map->numCols*map->map->numRows); + for (int i = 0; i < vectors->n; i++) { + vectors->data[i] = psVectorAllocEmpty (4, PS_TYPE_F32); + } + + for (int i = 0; i < x->n; i++) { + + int xRuff = x->data.F32[i] / map->binning->nXbin; + int yRuff = y->data.F32[i] / map->binning->nYbin; + + int bin = xRuff + yRuff*map->map->numCols; + + psVector *vector = vectors->data[bin]; + vector->data.F32[vector->n] = f->data.F32[i]; + psVectorExtend (vector, 4, 1); + } + + for (int iy = 0; iy < map->map->numRows; iy++) { + for (int ix = 0; ix < map->map->numCols; ix++) { + + // select the vector + psVector *vector = vectors->data[ix + iy*Nx]; + + // get the value + if (psVectorStats (stats, vector, NULL, NULL, 0)) { + mask->data.U8[iy][ix] = 0; + map->map->data.F32[iy][ix] = stats->robustMedian; // XXX select correct stats + } else { + mask->data.U8[iy][ix] = 1; + } + } + } + + psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff); + map->nGood = mask->numCols * mask->numRows - map->nBad - map->nPoor; + if (nBad > badFrac * mask->numCols * mask->numRows) { + return false: + } + + psImagePixelInterpolatePoor (map->map, mask, 0xff, state); + return true; +} Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h (revision 14721) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h (revision 14721) @@ -0,0 +1,41 @@ +/** @file psImageMap.c + * + * @brief Functions define a 2d coarse representation of a finer 2D field + * + * @ingroup Image + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-01 00:41:01 $ + * + * Copyright 2007 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_MAP_H +#define PS_IMAGE_MAP_H + +/// @addtogroup ImageOps Image Operations +/// @{ + +// a structure to describe the 2D variations of some quantity as a function of position the +// variation is represented as a psImage which covers the field also represented as a psImage. +// the map image pixels are superpixels of the field image. values in the field are determined +// by interpolating the map image. +typedef struct { + psStats *stats; + psImage *map; + psImage *field; + psImageMap *binning; + int nBad; + int nPoor; + int nGood; +} psImageMap; + +psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) PS_ATTR_MALLOC; + +// generate a psImageMap (or NULL) with the given number of superpixels in X and Y +psImageMap *psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac); + +/// @} +#endif // #ifndef PS_IMAGE_MAP_H Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (revision 14721) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (revision 14721) @@ -0,0 +1,147 @@ +/** @file psImagePixelInterpolate.c + * + * @brief Functions for interpolating bad pixels in images + * + * these functions test and set masked pixels in an image. These functions are complementary + * to the psImageInterpolate functions, which perform sub-pixel interpolation. Those + * functions require all pixels surrounding the sub-pixel interpolation to have values which + * are valid. These functions enable interpolation of complete missing pixels, potentially + * across large spans. + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-01 00:41:01 $ + * + * Copyright 2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include +#include +#include +#include + +#include "psAbort.h" +#include "psMemory.h" +#include "psError.h" +#include "psAssert.h" +#include "psString.h" +#include "psImage.h" +#include "psImageInterpolate.h" + +# define PS_IMAGE_ITER_STATE(XS,XE,YS,YE,N_MIN,TYPE) \ + nGood = 0; \ + for (jy = YS; jy <= YE; jy++) { \ + // stick to pixels in image grid \ + if (jy + iy < 0) { continue; } \ + if (jy + iy >= mask->numRows) { continue; } \ + for (jx = XS; jx <= XE; jx++) { \ + // stick to pixels in image grid \ + if (jx + ix < 0) { continue; } \ + if (jx + ix >= mask->numCols) { continue; } \ + // do not test self \ + if (!jx && jy) { continue; } \ + if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; } \ + nGood ++; \ + } \ + } \ + if (nGood >= N_MIN) { \ + nPoor ++; \ + result->data.S32[iy][ix] = TYPE; \ + continue; \ + } + +// count and mark pixels based on their potential for being interpolated. the input image is +// just the mask, the output image is the mask image in which the good, bad, and poor pixels +// are marked. +// XXX return in the image the state (2nd, which of 4 corners, etc)? +psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal) { + + psImage *result = psImageAlloc (mask->numCols, mask->numRows, PS_TYPE_S32); + psImageInit (result, 0); + + int nPoor = 0; + int nBad = 0; + + for (int iy = 0; iy < mask->numRows; iy++) { + for (int ix = 0; ix < mask->numCols; ix++) { + + // skip the good pixels (unmasked) + if (!(mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal)) { continue; } + + // examine the neighbors. If at least 6 of the 8 surrounding, or 3 of the 4 corner + // neighbors are valid, this is a pixel which can be interpolated. + + int nGood; + + // check for poor pixels + PS_IMAGE_ITER_STATE (-1,+1,-1,+1,6, PS_IMAGE_INTERPOLATE_CENTER); + PS_IMAGE_ITER_STATE (-1,+0,-1,+0,3, PS_IMAGE_INTERPOLATE_LL); + PS_IMAGE_ITER_STATE (-1,+0,+0,+1,3, PS_IMAGE_INTERPOLATE_UL); + PS_IMAGE_ITER_STATE (+0,+1,-1,+0,3, PS_IMAGE_INTERPOLATE_LR); + PS_IMAGE_ITER_STATE (+0,+1,+0,+1,3, PS_IMAGE_INTERPOLATE_UR); + + nBad ++; + result->data.S32[iy][ix] = PS_IMAGE_INTERPOLATE_BAD; + } + } + return result; +} + +// interpolate the poor pixels using the available options +bool psImagePixelInterpolatePoor (psImage *image, psImage *mask, psMaskType maskVal, psImage *state) { + + for (int iy = 0; iy < state->numRows; iy++) { + for (int ix = 0; ix < state->numCols; ix++) { + + switch (state->data.S32[iy][ix]) { + case PS_IMAGE_INTERPOLATE_GOOD: + // skip the good pixels + continue; + + case PS_IMAGE_INTERPOLATE_BAD: + // skip the bad pixels + continue; + + case PS_IMAGE_INTERPOLATE_CENTER: + // fit a quadratic to the valid neighbor pixels + // apply the quadratic to get the poor pixel value + psPolynomial2D *poly = psImageBicubeFit (image, mask, maskVal, ix, iy); + image->data.F32[iy][ix] = poly->coeff[0][0]; + continue; + + case PS_IMAGE_INTERPOLATE_LL: + // fit a plane to the valid neighbor pixels + // apply the plane to get the poor pixel value + psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 1, iy - 1); + image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[1][0] + poly->coeff[0][1]; + continue; + case PS_IMAGE_INTERPOLATE_LR: + // fit a plane to the valid neighbor pixels + // apply the plane to get the poor pixel value + psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 0, iy - 1); + image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[0][1]; + continue; + case PS_IMAGE_INTERPOLATE_UL: + // fit a plane to the valid neighbor pixels + // apply the plane to get the poor pixel value + psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 1, iy - 0); + image->data.F32[iy][ix] = poly->coeff[0][0] + poly->coeff[1][0]; + continue; + case PS_IMAGE_INTERPOLATE_UR: + // fit a plane to the valid neighbor pixels + // apply the plane to get the poor pixel value + psPolynomial2D *poly = psImageBilinearFit (image, mask, maskVal, ix - 0, iy - 0); + image->data.F32[iy][ix] = poly->coeff[0][0]; + continue; + } + } + } + return result; +} + Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h (revision 14721) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h (revision 14721) @@ -0,0 +1,37 @@ +/** @file psImagePixelInterpolate.c + * + * @brief Functions for interpolating bad pixels in images + * + * @ingroup Image + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-01 00:41:01 $ + * + * Copyright 2007 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_PIXEL_INTERPOLATE_H +#define PS_IMAGE_PIXEL_INTERPOLATE_H + +/// @addtogroup ImageOps Image Operations +/// @{ + +// XXX make these all bit values? +typedef enum { + PS_IMAGE_INTERPOLATE_GOOD = 0x00, + PS_IMAGE_INTERPOLATE_BAD = 0x01, + PS_IMAGE_INTERPOLATE_CENTER = 0x02, + PS_IMAGE_INTERPOLATE_CORNER = 0x04, + PS_IMAGE_INTERPOLATE_UR = 0x04, + PS_IMAGE_INTERPOLATE_UL = 0x05, + PS_IMAGE_INTERPOLATE_LR = 0x06, + PS_IMAGE_INTERPOLATE_LL = 0x07, +} psImagePixelInterpolateType; + +psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal); +bool psImagePixelInterpolatePoor (psImage *image, psImage *mask, psMaskType maskVal, psImage *state); + +/// @} +#endif // #ifndef PS_IMAGE_MAP_H Index: /branches/eam_branch_20070830/psLib/src/math/psPolynomialUtils.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/math/psPolynomialUtils.c (revision 14720) +++ /branches/eam_branch_20070830/psLib/src/math/psPolynomialUtils.c (revision 14721) @@ -135,5 +135,6 @@ } -// this function expect x,y in parent coords +// this function expects x,y in parent coords +// XXX add a mask, fit only the valid pixels psPolynomial2D *psImageBicubeFit(const psImage *image, int x, int y) { Index: /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am =================================================================== --- /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am (revision 14720) +++ /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am (revision 14721) @@ -22,4 +22,5 @@ tap_psImagePixelExtract \ tap_psImageInterpolate2 \ + tap_psImageMap \ tap_psImageMaskOps Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14721) +++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14721) @@ -0,0 +1,92 @@ +#include +#include +#include + +#include "tap.h" +#include "pstap.h" + +int SaveImage (psMetadata *header, psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + return (TRUE); +} + +int main (void) +{ + + // plan_tests(0); + + // *** tests to demonstrate the validity of the algorithm or concept *** + + // make a model for a well-sampled field of a simple function (f = ax + by + c) + { + // function is defined over the range 0-1000, 0-1000 + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 1000; + binning->nYfine = 1000; + binning->nXruff = 5; + binning->nYruff = 5; + + // generate a grid of test data points + psVector *x = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32); + + for (int ix = 0; ix < 1000; ix += 25) { + for (int iy = 0; iy < 1000; iy += 25) { + x->data.F32[x->n] = ix; + y->data.F32[y->n] = iy; + f->data.F32[f->n] = C00 + C10*ix + C01*iy; + psVectorExtend (x, 100, 1); + psVectorExtend (y, 100, 1); + psVectorExtend (f, 100, 1); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // allocate a map, but we have no field image to supply + bool status = psImageMapGenerate (map, x, y, f, 0.1); + psFree (binning); + + fprint (stderr, "nGood: %d\n", map->nGood); + fprint (stderr, "nPoor: %d\n", map->nPoor); + fprint (stderr, "nBad: %d\n", map->nBad); + + SaveImage (NULL, map->map, "map.fits"); + + psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + field->data.F32[iy][ix] = C00 + C10*ix + C01*iy; + } + } + SaveImage (NULL, field, "field.fits"); + + // measure difference between model (map) and data (field) + for (int ix = 0; ix < map->map->numCols; ix++) { + for (int iy = 0; iy < map->map->num; iy++) { + int xo = ix*map->binning->nXbin; + int yo = iy*map->binning->nYbin; + df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix]; + fprintf (stderr, "%d %d %f = %f - %f\n", ix, iy, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]); + } + } + } + + // make a model for a poorly-sampled field of a simple function (f = ax + by + c) + + // choose a model scale for a poorly-sampled field of a simple function (f = ax + by + c) + + // make a model for a well-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1) + + // make a model for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1) + + // choose a model scale for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1) + +}