Index: /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14725) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (revision 14726) @@ -7,6 +7,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 20:29:50 $ + * @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-03 20:29:07 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -63,15 +63,17 @@ psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK); + psImage *xCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32); + psImage *yCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32); // 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); + // we can do this by accumulating a vector of pixel indexes for each cell + psArray *pixelSets = psArrayAlloc (map->map->numCols*map->map->numRows); + for (int i = 0; i < pixelsSets->n; i++) { + pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32); } - + // associate each value with a cell for (int i = 0; i < x->n; i++) { - + // XXX use the psImageUnbin command to convert from x,y, to xRuff,yRuff? int xRuff = x->data.F32[i] / map->binning->nXbin; int yRuff = y->data.F32[i] / map->binning->nYbin; @@ -79,9 +81,10 @@ int bin = xRuff + yRuff*map->map->numCols; - psVector *vector = vectors->data[bin]; - vector->data.F32[vector->n] = f->data.F32[i]; + psVector *pixels = pixelSets->data[bin]; + pixels->data.S32[pixels->n] = i; psVectorExtend (vector, 4, 1); } + // accumulate the x,y coords for each point to calculate the mean position. int Nx = map->map->numCols; int Ny = map->map->numRows; @@ -89,6 +92,18 @@ for (int ix = 0; ix < Nx; ix++) { - // select the vector - psVector *vector = vectors->data[ix + iy*Nx]; + // select the vector. use psImageUnbin command to convert xRuff,yRuff to + // xFine,yFine + psVector *pixels = pixelSets->data[ix + iy*Nx]; + + psVector xCell = psVectorAlloc (pixels->n, PS_TYPE_F32); + psVector yCell = psVectorAlloc (pixels->n, PS_TYPE_F32); + psVector fCell = psVectorAlloc (pixels->n, PS_TYPE_F32); + + for (int i = 0; i < pixels->n; i++) { + bin = pixels->data.S32[i]; + xCell->data.F32[i] = x->data.F32[bin]; + yCell->data.F32[i] = y->data.F32[bin]; + fCell->data.F32[i] = f->data.F32[bin]; + } // reset the stats to avoid contamination from the previous loop @@ -96,13 +111,27 @@ // get the value - if (psVectorStats (map->stats, vector, NULL, NULL, 0)) { + // XXX need to supply a mask and skip the masked pixels when calculating the centroid + // this will not in general be properly weighted... + if (psVectorStats (map->stats, fCell, NULL, NULL, 0)) { mask->data.U8[iy][ix] = 0; // XXX ensure only one option is selected, or save both position and width map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options); + + // calculate the mean position and save: + psVectorStats (meanStat, xCell, NULL, NULL, 0); + xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); + psVectorStats (meanStat, yCell, NULL, NULL, 0); + yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); } else { mask->data.U8[iy][ix] = 1; } + + psFree (xCell); + psFree (yCell); + psFree (fCell); } } + + // at this point, for each map pixel, we have (f,x,y), or the pixel is masked. psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff); @@ -112,5 +141,12 @@ } + // fit the valid pixels to (0,1,2) order polynomials, interpolate values to the pixel center + // XXX I need to be careful about the pixel coordinates: center is 0,0 or 0.5, 0.5? + psImagePixelInterpolateCenters (map->map, xCoord, yCoord, state, mask, 0xff); + psFree (xCoord); + psFree (yCoord); + 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 14725) +++ /branches/eam_branch_20070830/psLib/src/imageops/psImagePixelInterpolate.c (revision 14726) @@ -11,6 +11,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 20:29:50 $ + * @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-03 20:29:07 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -163,4 +163,6 @@ break; } + // XXX should I use 1 1D polynomial fitting all unmasked pixels in the 3x3 grid? + // XXX that would automatically extend to regions where only 2 pixels are valid... case PS_IMAGE_INTERPOLATE_LL: { // fit a plane to the 3 pixels at (0,1),(1,0),(1,1), extend to pixel at (0,0) @@ -196,4 +198,123 @@ return true; } - +// interpolate the good pixels to their true centers +bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal) { + + assert (value->numCols == state->numCols); + assert (value->numRows == state->numRows); + assert (value->numCols == mask->numCols); + assert (value->numRows == mask->numRows); + + psImage *output = psImageAlloc (value->numCols, value->numRows, PS_TYPE_F32); + + // 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 *poly2D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); + poly2D->mask[2][2] = 1; + poly2D->mask[2][1] = 1; + poly2D->mask[1][2] = 1; + + // allocate a 2D polynomial to fit a plane to the valid neighbor pixels. + psPolynomial2D *poly1D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1); + poly2D->mask[1][1] = 1; + + 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_GOOD2: { + // XXX is there a fit-image-region function? + int n = 0; + for (int jy = -1; jy <= +1; jy++) { + // skip invalid pixels + if (jy + iy < 0) { continue; } + if (jy + iy >= value->numRows) { continue; } + for (int jx = -1; jx <= +1; jx++) { + // skip invalid pixels + if (jx + ix < 0) { continue; } + if (jx + ix >= value->numCols) { continue; } + // skip masked pixels + if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } + x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx]; + y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx]; + f->data.F32[n] = value->data.F32[iy+jy][ix+jx]; + // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx]; + n++; + } + } + // set vector lengths here + x->n = n; + y->n = n; + f->n = n; + // df->n = n; + // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y); + psVectorFitPolynomial2D (poly2D, NULL, 0xff, f, NULL, x, y); + // apply the fitted quadratic to get the poor pixel value + output->data.F32[iy][ix] = psPolynomial2DEval (poly2D, ix, iy); + break; } + + case PS_IMAGE_INTERPOLATE_GOOD1: { + // XXX is there a fit-image-region function? + int n = 0; + for (int jy = -1; jy <= +1; jy++) { + // skip invalid pixels + if (jy + iy < 0) { continue; } + if (jy + iy >= value->numRows) { continue; } + for (int jx = -1; jx <= +1; jx++) { + // skip invalid pixels + if (jx + ix < 0) { continue; } + if (jx + ix >= value->numCols) { continue; } + // skip masked pixels + if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; } + x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx]; + y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx]; + f->data.F32[n] = value->data.F32[iy+jy][ix+jx]; + // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx]; + n++; + } + } + // set vector lengths here + x->n = n; + y->n = n; + f->n = n; + // df->n = n; + // psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, df, x, y); + psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, NULL, x, y); + // apply the fitted quadratic to get the poor pixel value + output->data.F32[iy][ix] = psPolynomial2DEval (poly1D, ix, iy); + break; } + + case PS_IMAGE_INTERPOLATE_GOOD0: { + output->data.F32[iy][ix] = value->data.F32[iy+jy][ix+jx]; + break; } + + default: + psAbort("impossible case in __func__"); + } + } + } + + for (int iy = 0; iy < value->numRows; iy++) { + for (int ix = 0; ix < value->numCols; ix++) { + if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; } + value->data.F32[iy][ix] = output->data.F32[iy][ix]; + } + } + + psFree (x); + psFree (y); + psFree (f); + + psFree (poly2D); + psFree (poly1D); + + psFree (output); + return true; +} Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c =================================================================== --- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14725) +++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMap.c (revision 14726) @@ -92,4 +92,64 @@ // make a model for a well-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1) + { + // 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 += 10) { + for (int iy = 0; iy < 1000; iy += 10) { + x->data.F32[x->n] = ix; + y->data.F32[y->n] = iy; + f->data.F32[f->n] = C00 / (C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500) + 0.1); + 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 + // XXX this function needs to correct for the mean superpixel position + // which is sampled... + psImageMapGenerate (map, x, y, f, 0.1); + psFree (binning); + + 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"); + + 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*PS_SQR(ix-500) + C01*PS_SQR(iy-500) + 0.1); + } + } + 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->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]); + } + } + } // make a model for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)