Index: /trunk/psLib/test/imageops/.cvsignore =================================================================== --- /trunk/psLib/test/imageops/.cvsignore (revision 14925) +++ /trunk/psLib/test/imageops/.cvsignore (revision 14926) @@ -33,2 +33,6 @@ tap_psImageSmooth tap_psImageStructManip +tap_psImageMap + +tap_psImageMapFit +tap_psImageMapFit2 Index: /trunk/psLib/test/imageops/Makefile.am =================================================================== --- /trunk/psLib/test/imageops/Makefile.am (revision 14925) +++ /trunk/psLib/test/imageops/Makefile.am (revision 14926) @@ -22,4 +22,7 @@ tap_psImagePixelExtract \ tap_psImageInterpolate2 \ + tap_psImageMap \ + tap_psImageMapFit \ + tap_psImageMapFit2 \ tap_psImageMaskOps Index: /trunk/psLib/test/imageops/tap_psImageMap.c =================================================================== --- /trunk/psLib/test/imageops/tap_psImageMap.c (revision 14926) +++ /trunk/psLib/test/imageops/tap_psImageMap.c (revision 14926) @@ -0,0 +1,244 @@ +#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); +} + +# define C00 +0.5 +# define C01 +0.1 +# define C10 -0.2 + +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) + # if (0) + { + // 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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + 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 + // 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*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->numRows; iy++) { + int xo = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + 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]); + } + } + + psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel (ix + 0.5, iy + 0.5, map->map, map->binning); + } + } + SaveImage (NULL, model, "model.fits"); + } + # endif + + // 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 high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c) + # if (0) + { + // 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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + x->data.F32[x->n] = ix; + y->data.F32[y->n] = iy; + f->data.F32[f->n] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500); + 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] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500); + } + } + 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 = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + 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]); + } + } + } + # endif + + // make a model for a well-sampled field of a high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c) + # if (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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + x->data.F32[x->n] = ix; + y->data.F32[y->n] = iy; + f->data.F32[f->n] = C10*(ix-500) + PS_SQR(C01*(iy-500)); + 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] = C10*(ix-500) + PS_SQR(C01*(iy-500)); + } + } + 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 = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + 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]); + } + } + + + psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + // XXX the binned coordinates are probably off by 0.5 pix, or we are interpolating + // to the wrong binned coordinate. + // XXX fix edge cases + model->data.F32[iy][ix] = psImageUnbinPixel (ix + 0.5, iy + 0.5, map->map, map->binning); + } + } + SaveImage (NULL, model, "model.fits"); + } + # endif + + // 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) + +} Index: /trunk/psLib/test/imageops/tap_psImageMapFit.c =================================================================== --- /trunk/psLib/test/imageops/tap_psImageMapFit.c (revision 14926) +++ /trunk/psLib/test/imageops/tap_psImageMapFit.c (revision 14926) @@ -0,0 +1,843 @@ +#include +#include +#include + +#include "tap.h" +#include "pstap.h" + +// save function used to dump out test images while debugging algorithm +int SaveImage (psMetadata *header, psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + return (TRUE); +} + +# define TEST_4PT_0 1 +# define TEST_4PT_1 1 +# define TEST_4PT_2 1 +# define TEST_4PT_3 1 +# define TEST_6PT_0 1 +# define TEST_9PT_0 1 +# define TEST_9PT_1 1 +# define TEST_9PT_2 1 +# define TEST_9PT_3 1 +# define TEST_9PT_4 1 +# define TEST_9PT_5 1 + +int main (void) +{ + + plan_tests(219); + + // *** tests to demonstrate the validity of the algorithm or concept *** + + // test with unconstrained cell: 3x3 grid fitted to 8 points with simple slope and scale difference + # if (TEST_9PT_5) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 6; + binning->nYfine = 6; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (50, PS_TYPE_F32); + psVector *y = psVectorAlloc (50, PS_TYPE_F32); + psVector *f = psVectorAlloc (50, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (float ix = 1.0; ix < 6.0; ix += 2.0) { + for (float iy = 1.0; iy < 6.0; iy += 2.0) { + if ((ix == 1.0) && (iy == 1.0)) continue; + x->data.F32[n] = ix; + y->data.F32[n] = iy; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psImage *field = psImageAlloc(6, 6, PS_TYPE_F32); + for (int ix = 0; ix < 6; ix++) { + for (int iy = 0; iy < 6; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + if ((ix < 3) && (iy < 3)) { + is_float (model->data.F32[iy][ix], NAN, "model matches expected NaN"); + } else { + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs"); + } + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference + # if (TEST_9PT_4) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 6; + binning->nYfine = 6; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (50, PS_TYPE_F32); + psVector *y = psVectorAlloc (50, PS_TYPE_F32); + psVector *f = psVectorAlloc (50, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (float ix = 0.5; ix < 6.0; ix += 1.0) { + for (float iy = 0.5; iy < 6.0; iy += 1.0) { + x->data.F32[n] = ix; + y->data.F32[n] = iy; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psImage *field = psImageAlloc(6, 6, PS_TYPE_F32); + for (int ix = 0; ix < 6; ix++) { + for (int iy = 0; iy < 6; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference + # if (TEST_9PT_3) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 6; + binning->nYfine = 6; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (50, PS_TYPE_F32); + psVector *y = psVectorAlloc (50, PS_TYPE_F32); + psVector *f = psVectorAlloc (50, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (float ix = 1.0; ix < 6; ix += 2.0) { + for (float iy = 1.0; iy < 6; iy += 2.0) { + y->data.F32[n] = iy; + x->data.F32[n] = ix; + if ((ix == 1.0) && (iy == 1.0)) { + x->data.F32[n] = ix - 0.1; + } + if ((ix == 3.0) && (iy == 1.0)) { + y->data.F32[n] = iy - 0.1; + } + if ((ix == 5.0) && (iy == 3.0)) { + x->data.F32[n] = ix + 0.1; + } + if ((ix == 3.0) && (iy == 5.0)) { + y->data.F32[n] = iy + 0.1; + } + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psImage *field = psImageAlloc(6, 6, PS_TYPE_F32); + for (int ix = 0; ix < 6; ix++) { + for (int iy = 0; iy < 6; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference +# if (TEST_9PT_2) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 3; + binning->nYfine = 3; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (50, PS_TYPE_F32); + psVector *y = psVectorAlloc (50, PS_TYPE_F32); + psVector *f = psVectorAlloc (50, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (float ix = 0.0; ix < 3; ix += 1.0) { + for (float iy = 0.0; iy < 3; iy += 1.0) { + x->data.F32[n] = ix + 0.5; + y->data.F32[n] = iy + 0.5; + + if (ix == 0.0) { + x->data.F32[n] -= 0.1; + } + if (iy == 0.0) { + y->data.F32[n] -= 0.1; + } + if (ix == 2.0) { + x->data.F32[n] += 0.1; + } + if (iy == 2.0) { + y->data.F32[n] += 0.1; + } + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psImage *field = psImageAlloc(3, 3, PS_TYPE_F32); + for (int ix = 0; ix < 3; ix++) { + for (int iy = 0; iy < 3; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs"); + } + } + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } +# endif + + // still a simple test: 3x3 grid fitted to 9 points with simple slope and scale difference + # if (TEST_9PT_1) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 6; + binning->nYfine = 6; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (9, PS_TYPE_F32); + psVector *y = psVectorAlloc (9, PS_TYPE_F32); + psVector *f = psVectorAlloc (9, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (int ix = 1; ix < 6; ix += 2) { + for (int iy = 1; iy < 6; iy += 2) { + x->data.F32[n] = ix; + y->data.F32[n] = iy; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + + psImage *field = psImageAlloc(6, 6, PS_TYPE_F32); + for (int ix = 0; ix < 6; ix++) { + for (int iy = 0; iy < 6; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // very simple test: 3x3 grid fitted to 9 points with simple slope + # if (TEST_9PT_0) + { + // function is defined over the range 0-1000, 0-1000 + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 3; + binning->nYfine = 3; + binning->nXruff = 3; + binning->nYruff = 3; + + // generate a grid of test data points + psVector *x = psVectorAlloc (9, PS_TYPE_F32); + psVector *y = psVectorAlloc (9, PS_TYPE_F32); + psVector *f = psVectorAlloc (9, PS_TYPE_F32); + + int n = 0; + psImage *field = psImageAlloc(3, 3, PS_TYPE_F32); + for (int ix = 0; ix < 3; ix++) { + for (int iy = 0; iy < 3; iy++) { + x->data.F32[n] = ix + 0.5; + y->data.F32[n] = iy + 0.5; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + field->data.F32[iy][ix] = f->data.F32[n]; + n++; + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference + # if (TEST_6PT_0) + { + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 3; + binning->nYfine = 2; + binning->nXruff = 3; + binning->nYruff = 2; + + // generate a grid of test data points + psVector *x = psVectorAlloc (50, PS_TYPE_F32); + psVector *y = psVectorAlloc (50, PS_TYPE_F32); + psVector *f = psVectorAlloc (50, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + // place the measurement points exactly on the ruff reference pixel centers + int n = 0; + for (float ix = 0.5; ix < 3.0; ix += 1.0) { + for (float iy = 0.5; iy < 2.0; iy += 1.0) { + // x->data.F32[n] = ix; + // y->data.F32[n] = iy; + if ((ix == 0.5) && (iy == 0.5)) { + x->data.F32[n] = ix + 0.0; + y->data.F32[n] = iy + 0.1; // add in both points. + // f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + // n++; + // x->data.F32[n] = ix + 0.0; + // y->data.F32[n] = iy - 0.1; + } else { + x->data.F32[n] = ix; + y->data.F32[n] = iy; + } + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psImage *field = psImageAlloc(3, 2, PS_TYPE_F32); + for (int ix = 0; ix < 3; ix++) { + for (int iy = 0; iy < 2; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // more complex test: 2x2 grid fitted to 16 points with simple slope offset from grid centers + // this one uses points inset relative to the reference points so they are all visible + // to the reference + # if (TEST_4PT_3) + { + // function is defined over the range 0-1000, 0-1000 + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 4; + binning->nYfine = 4; + binning->nXruff = 2; + binning->nYruff = 2; + + // generate a grid of test data points + psVector *x = psVectorAlloc (16, PS_TYPE_F32); + psVector *y = psVectorAlloc (16, PS_TYPE_F32); + psVector *f = psVectorAlloc (16, PS_TYPE_F32); + + int n = 0; + // actual field is f = x + y, where x & y are the subpixel positions + for (float ix = 0.5; ix < 4; ix += 1.0) { + for (float iy = 0.5; iy < 4; iy += 1.0) { + x->data.F32[n] = ix; + y->data.F32[n] = iy; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + + psImage *field = psImageAlloc(4, 4, PS_TYPE_F32); + for (int ix = 0; ix < 4; ix++) { + for (int iy = 0; iy < 4; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 10*FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // more complex test: 2x2 grid fitted to 4 points with simple slope offset from grid centers + // this one uses points inset relative to the reference points so they are all visible + // to the reference + # if (TEST_4PT_2) + { + // function is defined over the range 0-1000, 0-1000 + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 2; + binning->nYfine = 2; + binning->nXruff = 2; + binning->nYruff = 2; + + // generate a grid of test data points + psVector *x = psVectorAlloc (5, PS_TYPE_F32); + psVector *y = psVectorAlloc (5, PS_TYPE_F32); + psVector *f = psVectorAlloc (5, PS_TYPE_F32); + + int n = 0; + psImage *field = psImageAlloc(2, 2, PS_TYPE_F32); + // actual field is f = x + y, where x & y are the subpixel positions + for (int ix = 0; ix < 2; ix++) { + for (int iy = 0; iy < 2; iy++) { + # if (1) + if (!ix && !iy) { + x->data.F32[n] = ix + 0.5 - 0.1; + y->data.F32[n] = iy + 0.5 - 0.0; +// f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; +// n++; +// x->data.F32[n] = ix + 0.5; +// y->data.F32[n] = iy + 0.5 - 0.1; + } else { + x->data.F32[n] = ix + 0.5; + y->data.F32[n] = iy + 0.5; + } + # endif + # if (0) + // offset points from centers + if (ix) { + x->data.F32[n] = ix + 0.5 - 0.1; + } else { + x->data.F32[n] = ix + 0.5 + 0.1; + } + if (iy) { + y->data.F32[n] = iy + 0.5 - 0.1; + } else { + y->data.F32[n] = iy + 0.5 + 0.1; + } + # endif + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + field->data.F32[iy][ix] = ix + 0.5 + iy + 0.5; + n++; + } + } + x->n = n; + y->n = n; + f->n = n; + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 5*FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // still a simple test: 2x2 grid fitted to 4 points with simple slope and scale difference + # if (TEST_4PT_1) + { + // function is defined over the range 0-1000, 0-1000 + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 4; + binning->nYfine = 4; + binning->nXruff = 2; + binning->nYruff = 2; + + // generate a grid of test data points + psVector *x = psVectorAlloc (4, PS_TYPE_F32); + psVector *y = psVectorAlloc (4, PS_TYPE_F32); + psVector *f = psVectorAlloc (4, PS_TYPE_F32); + + // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates + int n = 0; + for (int ix = 1; ix < 4; ix += 2) { + for (int iy = 1; iy < 4; iy += 2) { + x->data.F32[n] = ix; + y->data.F32[n] = iy; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + n++; + } + } + + psImage *field = psImageAlloc(4, 4, PS_TYPE_F32); + for (int ix = 0; ix < 4; ix++) { + for (int iy = 0; iy < 4; iy++) { + field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5); + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // simplest possible test: 2x2 grid fitted to 4 points with simple slope + # if (TEST_4PT_0) + { + // function is defined over the range 0-1000, 0-1000 + psMemId id = psMemGetId(); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 2; + binning->nYfine = 2; + binning->nXruff = 2; + binning->nYruff = 2; + + // generate a grid of test data points + psVector *x = psVectorAlloc (4, PS_TYPE_F32); + psVector *y = psVectorAlloc (4, PS_TYPE_F32); + psVector *f = psVectorAlloc (4, PS_TYPE_F32); + + int n = 0; + psImage *field = psImageAlloc(2, 2, PS_TYPE_F32); + for (int ix = 0; ix < 2; ix++) { + for (int iy = 0; iy < 2; iy++) { + x->data.F32[n] = ix + 0.5; + y->data.F32[n] = iy + 0.5; + f->data.F32[n] = x->data.F32[n] + y->data.F32[n]; + field->data.F32[iy][ix] = f->data.F32[n]; + n++; + } + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // scale defines both field and map image sizes (nXfine, nXruff) + psImageMap *map = psImageMapAlloc (NULL, binning, stats); + + // fit the data to the map + psImageMapFit (map, NULL, 0, x, y, f, NULL); + + psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32); + for (int ix = 0; ix < model->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs"); + } + } + + // SaveImage (NULL, map->map, "map.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, model, "model.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + return exit_status(); +} Index: /trunk/psLib/test/imageops/tap_psImageMapFit2.c =================================================================== --- /trunk/psLib/test/imageops/tap_psImageMapFit2.c (revision 14926) +++ /trunk/psLib/test/imageops/tap_psImageMapFit2.c (revision 14926) @@ -0,0 +1,257 @@ +#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); +} + +# define TEST1 0 +# define TEST2 1 +# define TEST3 0 + +# define C00 +0.0 +# define C01 -2.0 +# define C10 +1.0 + +int main (void) +{ + + plan_tests(1); + + // make a model for a well-sampled field of a simple function (f = ax + by + c) + # if (TEST1) + { + psMemId id = psMemGetId(); + + // function is defined over the range 0-1000, 0-1000 + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = 1000; + binning->nYfine = 1000; + binning->nXruff = 10; + binning->nYruff = 10; + + // 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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + 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); + + // fit the data to the map + psImageMapFit (map, x, y, f, NULL); + + 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; + } + } + + // 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 = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + is_float_tol (field->data.F32[yo][xo], map->map->data.F32[iy][ix], 1e-3, "model matches data"); + } + } + + // XXX test on the model or don't bother + psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + } + } + + // SaveImage (NULL, model, "model.fits"); + // SaveImage (NULL, field, "field.fits"); + // SaveImage (NULL, map->map, "map.fits"); + + psFree (model); + psFree (binning); + psFree (map); + psFree (stats); + psFree (field); + psFree (x); + psFree (y); + psFree (f); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + # endif + + // 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 high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c) + # if (TEST2) + { + // 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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + x->data.F32[x->n] = ix + 0,5; + y->data.F32[y->n] = iy + 0.5; + f->data.F32[f->n] = PS_SQR(x->data.F32[x->n]) + PS_SQR(y->data.F32[x->n]); + 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... + psImageMapFit (map, x, y, f, NULL); + psFree (binning); + + 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] = PS_SQR(ix+0.5) + PS_SQR(iy+0.5); + } + } + 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 = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + 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]); + } + } + + // XXX test on the model or don't bother + psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + } + } + SaveImage (NULL, model, "model.fits"); + } + # endif + + // make a model for a well-sampled field of a high-order polynomial function (f = (a(x-xo)^2 + b(y-yo)^2) + c) + # if (TEST3) + { + // 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 += 20) { + for (int iy = 0; iy < 1000; iy += 20) { + x->data.F32[x->n] = ix; + y->data.F32[y->n] = iy; + f->data.F32[f->n] = C10*(ix-500) + PS_SQR(C01*(iy-500)); + 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] = C10*(ix-500) + PS_SQR(C01*(iy-500)); + } + } + 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 = psImageBinningGetFineX(map->binning, ix + 0.5); + int yo = psImageBinningGetFineY(map->binning, iy + 0.5); + 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]); + } + } + + + psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32); + for (int ix = 0; ix < 1000; ix++) { + for (int iy = 0; iy < 1000; iy++) { + // XXX the binned coordinates are probably off by 0.5 pix, or we are interpolating + // to the wrong binned coordinate. + // XXX fix edge cases + model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning); + } + } + SaveImage (NULL, model, "model.fits"); + } + # endif + + // 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) + + return exit_status(); +} Index: /trunk/psLib/test/imageops/tap_psImageStructManip.c =================================================================== --- /trunk/psLib/test/imageops/tap_psImageStructManip.c (revision 14925) +++ /trunk/psLib/test/imageops/tap_psImageStructManip.c (revision 14926) @@ -10,6 +10,6 @@ * psImageTrim() * -* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-06-04 20:25:32 $ +* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-09-20 23:56:10 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -33,166 +33,188 @@ psRegion region2 = psRegionSet(numCols/4, numCols/4+numCols/2, numRows/4, numRows/4+numRows/2); - psImage* original = psImageAlloc(numCols,numRows,PS_TYPE_U32); - for (psS32 row=0;rowdata.F32[row][col] = row*1000+col; - } - } - - memcpy(&preSubsetStruct, original, sizeof(psImage)); - psImage* subset2 = psImageSubset(original, region2); - ok(subset2, "psImageSubset() returned non-NULL (subset2)"); - skip_start(subset2 == NULL, 25, "Skipping tests because psImageSubset() returned NULL"); - psImage* subset3 = psImageSubset(original, region1); - ok(subset3, "psImageSubset() returned non-NULL (subset3)"); - skip_start(subset3 == NULL, 24, "Skipping tests because psImageSubset() returned NULL"); - - // Verify the returned psImage structure members nrow and ncol are equal to - // the input parameter nrow and ncol respectively - ok(subset2->numCols == numCols/2 && subset2->numRows == numRows/2, - "psImageSubset output size set properly"); - ok(subset3->numCols == numCols/2 && subset3->numRows == numRows/2, - "psImageSubset output size set properly"); - - // Verify the returned psImage structure contains expected values in the - // row member, if the input psImage structure image contains known values - bool errorFlag = false; - for (psS32 row=0;rowdata.U32[row][col] != original->data.U32[row+numRows/4][col+numCols/4]) { - diag("psImageSubset output #1 was wrong at %dx%d (%d vs %d).", - row,col,subset2->data.U32[row][col], original->data.U32[row+numRows/4][col+numCols/4]); - errorFlag = true; - } - if (subset3->data.U32[row][col] != original->data.U32[row][col]) { - diag("psImageSubset output #1 was wrong at %dx%d (%d vs %d).", - row,col,subset2->data.U32[row][col], original->data.U32[row][col]); - errorFlag = true; - } - } - } - ok(!errorFlag, "psImageSubset() produced the expected values"); - - // Verify the returned psImage structure member type is equal to the input - // psImage structure member type - ok(subset2->type.type == PS_TYPE_U32, "psImageSubset() type was correct (subset2)"); - ok(subset3->type.type == PS_TYPE_U32, "psImageSubset() type was correct (subset3)"); - - // Verify the returned psImage structure members row0 and col0 are equal to - // the input parameters row0 and col0 respectively - ok(subset2->col0 == numCols/4 && subset2->row0 == numRows/4, - "psImageSubset() set col0/row0 correctly (subset2)"); - ok(subset3->col0 == 0 && subset3->row0 == 0, - "psImageSubset() set col0/row0 correctly (subset3)"); - - // Verify the returned psImage structure member parent is equal to the - // input psImage structure pointer image - ok(subset2->parent == original && subset3->parent == original, "psImageSubset() set ->parent correctly"); - - // Verify the returned psImage structure member children is null - ok(subset2->children == NULL && subset3->children == NULL, "psImageSubset() set ->children correctly"); - - // Verify the input psImage structure image only has the following members - // changed: 1) Nchildren is increased by one. 2) parent contains pointer psImage structure - // out at parent[Nchildren-1] - ok(original->children != NULL && original->children->n == 2, "psImageSubset did increment number of children by one per subset."); - ok(original->children->data[0] == subset2 && original->children->data[1] == subset3, - "psImageSubset did properly store the children pointers."); - - // Verify the returned psImage structure pointer is null and program - // execution doesn't stop, if the input parameter image is null. - // Also verified the input psImage structure is not modified - // An error should follow... - // XXX: Verify error - psImage* subset1 = psImageSubset(NULL,region1); - ok(subset1 == NULL, "psImageSubset returned NULL when input image was NULL."); - - // Verify the returned psImage structure pointer is null and program - // execution doesn't stop, if the input parameters nrow and/or ncol are zero. - // Also verify input psImage structure is not modified - // An error should follow... - // XXX: Verify error - { - memcpy(&preSubsetStruct,original,sizeof(psImage)); - subset1 = psImageSubset(original, psRegionSet(0,numCols/2,numRows/2,numRows/2)); - ok(subset1 == NULL, "psImageSubset returned NULL when numRows=0."); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original,psRegionSet(numCols/2,numCols/2,0,numRows/2)); - ok(subset1 == NULL, "psImageSubset returned NULL when numCols=0."); - ok(memcmp(original,&preSubsetStruct,sizeof(psImage)) == 0, - "psImageSubset didn't change the original struct though it failed to subset."); - } - - // Verify the returned psImage structure pointer is null and program - // execution doesn't stop, if the input parameters row0 and col0 are not within - // the range of values of psImage structure image - // An error should follow... - // XXX: Verify error - { - subset1 = psImageSubset(original, psRegionSet(0,numCols/2, 0,numRows*2)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset origin was outside of image (via cols)"); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original,psRegionSet(0,numCols*2,0,numRows/2)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset origin was outside of image (via rows)"); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original, psRegionSet(-1,numCols/2,0,numRows/2)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset origin was outside of image (col0=-1)"); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original, psRegionSet(0,numCols/2,-1,numRows/2)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset origin was outside of image (row0=-1)"); - } - - // Verify the returned psImage structure pointer is null and program - // execution doesn't stop if the input parameters nrow, ncol, row0 and col0 - // specify a range of data not within the input psImage structure image. Also - // verify the input psImage structure is not modified - // An error should follow... - // XXX: Verify error - { - subset1 = psImageSubset(original,psRegionSet(0,numCols/2,0,numRows+1)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset was outside of image (via rows)"); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original, psRegionSet(0,numCols+1,0,numRows/2)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset was outside of image (via cols)"); - // An error should follow... - // XXX: Verify error - subset1 = psImageSubset(original,psRegionSet(0,numCols+1,0,numRows+1)); - ok(subset1 == NULL, - "psImageSubset returned NULL when subset was outside of image (via row+cols)"); - } - - // psImageFreeChildren shall deallocate any children images of a - // psImage structure - memcpy(&preSubsetStruct,original,sizeof(psImage)); - psImageFreeChildren(original); - - // Verify the returned psImage structure member Nchildren is set to zero. - // XXX: This doesn't make sense. - ok(original->children == NULL || original->children->n == 0, - "psImageFreeChildren didn't set children to NULL"); - - //Verify the returned psImage structure members type, nrow, ncol, row0, col0, rows - // and parent are not modified. - ok(preSubsetStruct.numRows == original->numRows && - preSubsetStruct.numCols == original->numCols && - preSubsetStruct.row0 == original->row0 && - preSubsetStruct.col0 == original->col0, - "psImageFreeChildren modified parent's non-children elements"); - - skip_end(); - skip_end(); - psFree(original); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + // test basic subset creation + { + psImage* original = psImageAlloc(numCols,numRows,PS_TYPE_U32); + for (psS32 row=0;rowdata.F32[row][col] = row*1000+col; + } + } + + // XXX this is not being used in this section + memcpy(&preSubsetStruct, original, sizeof(psImage)); + psImage* subset1 = psImageSubset(original, region2); + ok(subset1, "psImageSubset() returned non-NULL (subset1)"); + skip_start(subset1 == NULL, 25, "Skipping tests because psImageSubset() returned NULL"); + psImage* subset2 = psImageSubset(original, region1); + ok(subset2, "psImageSubset() returned non-NULL (subset2)"); + skip_start(subset2 == NULL, 24, "Skipping tests because psImageSubset() returned NULL"); + + // Verify the returned psImage structure members nrow and ncol are equal to + // the input parameter nrow and ncol respectively + ok(subset1->numCols == numCols/2 && subset1->numRows == numRows/2, + "psImageSubset output size set properly"); + ok(subset2->numCols == numCols/2 && subset2->numRows == numRows/2, + "psImageSubset output size set properly"); + + // Verify the returned psImage structure contains expected values in the + // row member, if the input psImage structure image contains known values + bool errorFlag = false; + for (psS32 row=0;rowdata.U32[row][col] != original->data.U32[row+numRows/4][col+numCols/4]) { + diag("psImageSubset output #1 was wrong at %dx%d (%d vs %d).", + row,col,subset1->data.U32[row][col], original->data.U32[row+numRows/4][col+numCols/4]); + errorFlag = true; + } + if (subset2->data.U32[row][col] != original->data.U32[row][col]) { + diag("psImageSubset output #1 was wrong at %dx%d (%d vs %d).", + row,col,subset1->data.U32[row][col], original->data.U32[row][col]); + errorFlag = true; + } + } + } + ok(!errorFlag, "psImageSubset() produced the expected values"); + + // Verify the returned psImage structure member type is equal to the input + // psImage structure member type + ok(subset1->type.type == PS_TYPE_U32, "psImageSubset() type was correct (subset1)"); + ok(subset2->type.type == PS_TYPE_U32, "psImageSubset() type was correct (subset2)"); + + // Verify the returned psImage structure members row0 and col0 are equal to + // the input parameters row0 and col0 respectively + ok(subset1->col0 == numCols/4 && subset1->row0 == numRows/4, + "psImageSubset() set col0/row0 correctly (subset1)"); + ok(subset2->col0 == 0 && subset2->row0 == 0, + "psImageSubset() set col0/row0 correctly (subset2)"); + + // Verify the returned psImage structure member parent is equal to the + // input psImage structure pointer image + ok(subset1->parent == original && subset2->parent == original, "psImageSubset() set ->parent correctly"); + + // Verify the returned psImage structure member children is null + ok(subset1->children == NULL && subset2->children == NULL, "psImageSubset() set ->children correctly"); + + // Verify the input psImage structure image only has the following members + // changed: 1) Nchildren is increased by one. 2) parent contains pointer psImage structure + // out at parent[Nchildren-1] + ok(original->children != NULL && original->children->n == 2, "psImageSubset did increment number of children by one per subset."); + ok(original->children->data[0] == subset1 && original->children->data[1] == subset2, + "psImageSubset did properly store the children pointers."); + + psFree (subset1); + psFree (subset2); + skip_end(); + skip_end(); + psFree (original); + } + + // test free in reverse order from above + { + psImage* original = psImageAlloc(numCols,numRows,PS_TYPE_U32); + for (psS32 row=0;rowdata.F32[row][col] = row*1000+col; + } + } + + // XXX this is not being used in this section + memcpy(&preSubsetStruct, original, sizeof(psImage)); + psImage* subset1 = psImageSubset(original, region2); + ok(subset1, "psImageSubset() returned non-NULL (subset1)"); + psImage* subset2 = psImageSubset(original, region1); + ok(subset2, "psImageSubset() returned non-NULL (subset2)"); + + psFree (original); + psFree (subset1); + psFree (subset2); + } + + { + psImage* original = psImageAlloc(numCols,numRows,PS_TYPE_U32); + for (psS32 row=0;rowdata.F32[row][col] = row*1000+col; + } + } + + // XXX this is not being used in this section + memcpy(&preSubsetStruct, original, sizeof(psImage)); + + // Verify the returned psImage structure pointer is null and program + // execution doesn't stop, if the input parameter image is null. + // Also verified the input psImage structure is not modified + // An error should follow... + // XXX: Verify error + psImage* subset1 = psImageSubset(NULL,region1); + ok(subset1 == NULL, "psImageSubset returned NULL when input image was NULL."); + + // Verify the returned psImage structure pointer is null and program + // execution doesn't stop, if the input parameters nrow and/or ncol are zero. + // Also verify input psImage structure is not modified + // An error should follow... + // XXX: Verify error + { + memcpy(&preSubsetStruct,original,sizeof(psImage)); + subset1 = psImageSubset(original, psRegionSet(0,numCols/2,numRows/2,numRows/2)); + ok(subset1 == NULL, "psImageSubset returned NULL when numRows=0."); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original,psRegionSet(numCols/2,numCols/2,0,numRows/2)); + ok(subset1 == NULL, "psImageSubset returned NULL when numCols=0."); + ok(memcmp(original,&preSubsetStruct,sizeof(psImage)) == 0, + "psImageSubset didn't change the original struct though it failed to subset."); + } + + // Verify the returned psImage structure pointer is null and program + // execution doesn't stop, if the input parameters row0 and col0 are not within + // the range of values of psImage structure image + // An error should follow... + // XXX: Verify error + { + subset1 = psImageSubset(original, psRegionSet(0,numCols/2, 0,numRows*2)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset origin was outside of image (via cols)"); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original,psRegionSet(0,numCols*2,0,numRows/2)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset origin was outside of image (via rows)"); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original, psRegionSet(-1,numCols/2,0,numRows/2)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset origin was outside of image (col0=-1)"); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original, psRegionSet(0,numCols/2,-1,numRows/2)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset origin was outside of image (row0=-1)"); + } + + // Verify the returned psImage structure pointer is null and program + // execution doesn't stop if the input parameters nrow, ncol, row0 and col0 + // specify a range of data not within the input psImage structure image. Also + // verify the input psImage structure is not modified + // An error should follow... + // XXX: Verify error + { + subset1 = psImageSubset(original,psRegionSet(0,numCols/2,0,numRows+1)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset was outside of image (via rows)"); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original, psRegionSet(0,numCols+1,0,numRows/2)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset was outside of image (via cols)"); + // An error should follow... + // XXX: Verify error + subset1 = psImageSubset(original,psRegionSet(0,numCols+1,0,numRows+1)); + ok(subset1 == NULL, + "psImageSubset returned NULL when subset was outside of image (via row+cols)"); + } + psFree(original); + } + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); } }