IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14864


Ignore:
Timestamp:
Sep 17, 2007, 10:17:21 AM (19 years ago)
Author:
magnier
Message:

added errors to psImageMap, NaN for unconstrained pixels in output map

Location:
branches/eam_branch_20070830/psLib
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c

    r14774 r14864  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-07 20:19:35 $
     9 *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-17 20:17:21 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    3939
    4040    psFree (map->map);
    41     // psFree (map->field); XXX ??? this should be freed here, but that causes an error...
     41    psFree (map->error);
     42    psFree (map->field); // XXX ??? this should be freed here, but that causes an error...
    4243    psFree (map->stats);
    4344    psFree (map->binning);
    … …  
    5758    map->field   = psMemIncrRefCounter (field);
    5859    map->stats   = psMemIncrRefCounter (stats);
     60
    5961    map->map     = psImageAlloc (binning->nXruff, binning->nYruff, PS_TYPE_F32);
     62    psImageInit (map->map, 0.0);
     63
     64    map->error   = psImageAlloc (binning->nXruff, binning->nYruff, PS_TYPE_F32);
     65    psImageInit (map->error, 0.0);
    6066
    6167    psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
    … …  
    7379
    7480    psImageRecycle (map->map, nXruff, nYruff, PS_TYPE_F32);
     81    psImageRecycle (map->error, nXruff, nYruff, PS_TYPE_F32);
    7582
    7683    psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
    … …  
    257264}
    258265
     266// x,y are in fractional pixel coords of the fine image (pixel center: 0.5)
    259267double psImageMapEval (psImageMap *map, float x, float y) {
    260268
    261269    double result;
    262270
    263     result = psImageUnbinPixel(x, y, map->map, map->binning);
     271    result = psImageUnbinPixel_V2(x, y, map->map, map->binning);
    264272
    265273    return result;
    … …  
    276284
    277285    for (int i = 0; i < x->n; i++) {
    278         result->data.F32[i] = psImageUnbinPixel(x->data.F32[i], y->data.F32[i], map->map, map->binning);
     286        result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
    279287    }
    280288

  • branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h

    r14859 r14864  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-15 19:43:43 $
     9 *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-17 20:17:21 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    2626    psStats *stats;
    2727    psImage *map;
     28    psImage *error;
    2829    psImage *field;
    2930    psImageBinning *binning;

  • branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c

    r14862 r14864  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-17 01:14:52 $
     9 *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-17 20:17:21 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    2222
    2323// XXX for testing
    24 // #include "psFits.h"
    25 // #include "psFitsImage.h"
     24#include "psFits.h"
     25#include "psFitsImage.h"
    2626
    2727#include "psAssert.h"
    … …  
    238238    // fprintf (stderr, "Total: %f\n", Total);
    239239
    240     // write 1s in unused diagonal elements
    241     // XXX for test only
    242     # if (0)
    243     for (int n = 0; n < Nx; n+=Nx-1) {
    244         for (int m = 0; m < Ny; m++) {
    245             I = n + Nx * m;
    246             A->data.F32[I][I] = 1.0;
    247         }
    248     }
    249     for (int n = 0; n < Nx; n++) {
    250         for (int m = 0; m < Ny; m+=Ny-1) {
    251             I = n + Nx * m;
    252             A->data.F32[I][I] = 1.0;
    253         }
    254     }
    255     # endif
    256 
    257     // psFits *fits = psFitsOpen ("Agj.fits", "w");
    258     // psFitsWriteImage (fits, NULL, A, 0, NULL);
    259     // psFitsClose (fits);
     240    // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
     241    psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8);
     242    psVectorInit (Empty, 0);
     243    for (int i = 0; i < Nx*Ny; i++) {
     244        if (A->data.F32[i][i] == 0.0) {
     245            Empty->data.S8[i] = 1;
     246            for (int j = 0; j < Nx*Ny; j++) {
     247                A->data.F32[i][j] = 0.0;
     248                A->data.F32[j][i] = 0.0;
     249            }
     250            A->data.F32[i][i] = 1.0;
     251            B->data.F32[i] = 0.0;
     252        }
     253    }
     254
     255    psFits *fits = psFitsOpen ("Agj.fits", "w");
     256    psFitsWriteImage (fits, NULL, A, 0, NULL);
     257    psFitsClose (fits);
     258
     259    psImage *vector = psImageAlloc (1, B->n, PS_TYPE_F32);
     260    for (int n = 0; n < B->n; n++) {
     261        vector->data.F32[0][n] = B->data.F32[n];
     262    }
     263
     264    fits = psFitsOpen ("Bgj.fits", "w");
     265    psFitsWriteImage (fits, NULL, vector, 0, NULL);
     266    psFitsClose (fits);
     267    psFree (vector);
    260268
    261269    if (!psMatrixGJSolveF32(A, B)) {
    … …  
    267275    }
    268276   
     277    // set bad values to NaN
     278    for (int i = 0; i < Nx*Ny; i++) {
     279        if (Empty->data.S8[i]) {
     280            B->data.F32[i] = NAN;
     281        }
     282    }
     283
     284
    269285    for (int n = 0; n < Nx; n++) {
    270286        for (int m = 0; m < Ny; m++) {
    271287            I = n + Nx * m;
    272288            map->map->data.F32[m][n] = B->data.F32[I];
     289            map->error->data.F32[m][n] = sqrt(A->data.F32[I][I]);
    273290        }
    274291    }
    … …  
    276293    psFree (A);
    277294    psFree (B);
     295    psFree (Empty);
    278296
    279297    return true;

  • branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c

    r14863 r14864  
    2525# define TEST_9PT_3 1
    2626# define TEST_9PT_4 1
     27# define TEST_9PT_5 1
    2728
    2829int main (void)
    2930{
    3031
    31     plan_tests(1);
     32    plan_tests(219);
    3233   
    3334    // *** tests to demonstrate the validity of the algorithm or concept ***
     35
     36    // test with unconstrained cell: 3x3 grid fitted to 8 points with simple slope and scale difference
     37    # if (TEST_9PT_5)
     38    {
     39        psMemId id = psMemGetId();
     40
     41        psImageBinning *binning = psImageBinningAlloc();
     42        binning->nXfine = 6;
     43        binning->nYfine = 6;
     44        binning->nXruff = 3;
     45        binning->nYruff = 3;
     46
     47        // generate a grid of test data points
     48        psVector *x = psVectorAlloc (50, PS_TYPE_F32);
     49        psVector *y = psVectorAlloc (50, PS_TYPE_F32);
     50        psVector *f = psVectorAlloc (50, PS_TYPE_F32);
     51
     52        // the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
     53        // place the measurement points exactly on the ruff reference pixel centers
     54        int n = 0;
     55        for (float ix = 1.0; ix < 6.0; ix += 2.0) {
     56            for (float iy = 1.0; iy < 6.0; iy += 2.0) {
     57                if ((ix == 1.0) && (iy == 1.0)) continue;
     58                x->data.F32[n] = ix;
     59                y->data.F32[n] = iy;
     60                f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
     61                n++;
     62            }
     63        }
     64        x->n = n;
     65        y->n = n;
     66        f->n = n;
     67
     68        psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
     69        for (int ix = 0; ix < 6; ix++) {
     70            for (int iy = 0; iy < 6; iy++) {
     71                field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
     72            }
     73        }
     74
     75        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     76
     77        // scale defines both field and map image sizes (nXfine, nXruff)
     78        psImageMap *map = psImageMapAlloc (NULL, binning, stats);
     79
     80        // fit the data to the map
     81        psImageMapFit (map, x, y, f, NULL);
     82
     83        psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
     84        for (int ix = 0; ix < model->numCols; ix++) {
     85            for (int iy = 0; iy < model->numRows; iy++) {
     86                model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
     87                if ((ix < 3) && (iy < 3)) {
     88                    is_float (model->data.F32[iy][ix], NAN, "model matches expected NaN");
     89                } else {
     90                    is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs");
     91                }
     92            }
     93        }
     94
     95        SaveImage (NULL, map->map, "map.fits");
     96        SaveImage (NULL, field, "field.fits");
     97        SaveImage (NULL, model, "model.fits");
     98
     99        psFree (model);
     100        psFree (binning);
     101        psFree (map);
     102        psFree (stats);
     103        psFree (field);
     104        psFree (x);
     105        psFree (y);
     106        psFree (f);
     107        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     108    }
     109    # endif
    34110
    35111    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
Note: See TracChangeset for help on using the changeset viewer.