Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14724) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14725) @@ -7,6 +7,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 02:03:58 $ + * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-02 20:29:50 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -112,5 +112,5 @@ } - psImagePixelInterpolatePoor (map->map, mask, 0xff, state); + psImagePixelInterpolatePoor (map->map, state, mask, 0xff); return true; } Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (revision 14724) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (revision 14725) @@ -11,6 +11,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 02:42:55 $ + * @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-02 20:29:50 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -62,7 +62,6 @@ // 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)? +// just the mask, the output image contains enum values which define the type of interpolation which +// can be performed psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal) { @@ -99,5 +98,23 @@ // interpolate the poor pixels using the available options -bool psImagePixelInterpolatePoor (psImage *image, psImage *mask, psMaskType maskVal, psImage *state) { +bool psImagePixelInterpolatePoor (psImage *image, psImage *state, psImage *mask, psMaskType maskVal) { + + assert (image->numCols == state->numCols); + assert (image->numRows == state->numRows); + assert (image->numCols == mask->numCols); + assert (image->numRows == mask->numRows); + + // allocate the vectors for the 2nd order fit below + psVector *f = psVectorAlloc (9, PS_TYPE_F32); + // XXX if we add the weight above, include df + // psVector *df = psVectorAlloc (9, PS_TYPE_F32); + psVector *x = psVectorAlloc (9, PS_TYPE_F32); + psVector *y = psVectorAlloc (9, PS_TYPE_F32); + + // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels. + psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); + poly->mask[2][2] = 1; + poly->mask[2][1] = 1; + poly->mask[1][2] = 1; for (int iy = 0; iy < state->numRows; iy++) { @@ -114,41 +131,33 @@ case PS_IMAGE_INTERPOLATE_CENTER: { - // fit a quadratic to the valid neighbor pixels - // apply the quadratic to get the poor pixel value - psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly->mask[2][2] = 1; - poly->mask[2][1] = 1; - poly->mask[1][2] = 1; - psVector *f = psVectorAlloc (9, PS_TYPE_F32); - psVector *df = psVectorAlloc (9, PS_TYPE_F32); - psVector *x = psVectorAlloc (9, PS_TYPE_F32); - psVector *y = psVectorAlloc (9, PS_TYPE_F32); - + // XXX is there a fit-image-region function? int n = 0; for (int jy = -1; jy <= +1; jy++) { - // skip invalid pixels (jx < 0, >= Nx), etc + // skip invalid pixels if (jy + iy < 0) { continue; } if (jy + iy >= image->numRows) { continue; } for (int jx = -1; jx <= +1; jx++) { - // skip invalid pixels (jx < 0, >= Nx), etc + // skip invalid pixels if (jx + ix < 0) { continue; } - if (jx + ix >= mask->numCols) { continue; } - /* skip self */ + if (jx + ix >= image->numCols) { continue; } + // skip self if (!jx && !jy) { continue; } // skip masked pixels - if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } \ - f->data.F32[n] = image->data.F32[iy+jy][ix+jx]; - // df->data.F32[n] = weight->data.F32[jy][jx]; + if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } x->data.F32[n] = jx; y->data.F32[n] = jy; + f->data.F32[n] = image->data.F32[iy+jy][ix+jx]; + // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx]; n++; } } - // XXX set vector lengths here + // set vector lengths here x->n = n; y->n = n; f->n = n; // df->n = n; - psVectorFitPolynomial2D (poly, mask, maskValue, f, df, x, y); + // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y); + psVectorFitPolynomial2D (poly, NULL, 0xff, f, NULL, x, y); + // apply the fitted quadratic to get the poor pixel value image->data.F32[iy][ix] = poly->coeff[0][0]; break; } @@ -179,4 +188,10 @@ } } + + psFree (x); + psFree (y); + psFree (f); + + psFree (poly); return true; } Index: /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h (revision 14724) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.h (revision 14725) @@ -7,6 +7,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-01 00:41:01 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-02 20:29:50 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -32,5 +32,5 @@ psImage *psImagePixelInterpolateState (int *nBad, int *nPoor, psImage *mask, psMaskType maskVal); -bool psImagePixelInterpolatePoor (psImage *image, psImage *mask, psMaskType maskVal, psImage *state); +bool psImagePixelInterpolatePoor (psImage *image, psImage *state, psImage *mask, psMaskType maskVal); /// @} Index: /branches/eam_branch_20070830/psLib/src/pslib_strict.h =================================================================== --- /branches/eam_branch_20070830/psLib/src/pslib_strict.h (revision 14724) +++ /branches/eam_branch_20070830/psLib/src/pslib_strict.h (revision 14725) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-04-04 22:42:02 $ +* @version $Revision: 1.29.8.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-09-02 20:29:50 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -49,4 +49,5 @@ #include "psImagePixelExtract.h" #include "psImagePixelManip.h" +#include "psImagePixelInterpolate.h" #include "psImageStats.h" #include "psImageStructManip.h" @@ -54,4 +55,5 @@ #include "psImageBinning.h" #include "psImageUnbin.h" +#include "psImageMap.h" #include "psImageJpeg.h" Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14724) +++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14725) @@ -13,4 +13,8 @@ return (TRUE); } + +# define C00 +0.5 +# define C01 +0.1 +# define C10 -0.2 int main (void) @@ -35,6 +39,6 @@ psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32); - for (int ix = 0; ix < 1000; ix += 25) { - for (int iy = 0; iy < 1000; iy += 25) { + for (int ix = 0; ix < 1000; ix += 10) { + for (int iy = 0; iy < 1000; iy += 10) { x->data.F32[x->n] = ix; y->data.F32[y->n] = iy; @@ -52,10 +56,12 @@ // allocate a map, but we have no field image to supply - bool status = psImageMapGenerate (map, x, y, f, 0.1); + // XXX this function needs to correct for the mean superpixel position + // which is sampled... + 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); + fprintf (stderr, "nGood: %d\n", map->nGood); + fprintf (stderr, "nPoor: %d\n", map->nPoor); + fprintf (stderr, "nBad: %d\n", map->nBad); SaveImage (NULL, map->map, "map.fits"); @@ -71,9 +77,10 @@ // 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]); + for (int iy = 0; iy < map->map->numRows; iy++) { + + int xo = (ix + 0.5)*map->binning->nXbin; + int yo = (iy + 0.5)*map->binning->nYbin; + float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix]; + fprintf (stderr, "%d %d -> %d %d : %f = %f - %f\n", ix, iy, xo, yo, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]); } }