IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14726


Ignore:
Timestamp:
Sep 3, 2007, 10:29:07 AM (19 years ago)
Author:
magnier
Message:

updates to image map code, tests

Location:
branches/eam_branch_20070830/psLib
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r14725 r14726  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-02 20:29:50 $
     9 *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-03 20:29:07 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    6363
    6464    psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK);
     65    psImage *xCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32);
     66    psImage *yCoord = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_F32);
    6567
    6668    // accumulate the values for each map pixel
    6769   
    68     // we can do this by accumulating a sub-vector for each cell
    69     psArray *vectors = psArrayAlloc (map->map->numCols*map->map->numRows);
    70     for (int i = 0; i < vectors->n; i++) {
    71         vectors->data[i] = psVectorAllocEmpty (4, PS_TYPE_F32);
     70    // we can do this by accumulating a vector of pixel indexes for each cell
     71    psArray *pixelSets = psArrayAlloc (map->map->numCols*map->map->numRows);
     72    for (int i = 0; i < pixelsSets->n; i++) {
     73        pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
    7274    }
    73 
     75    // associate each value with a cell
    7476    for (int i = 0; i < x->n; i++) {
    75        
     77        // XXX use the psImageUnbin command to convert from x,y, to xRuff,yRuff?
    7678        int xRuff = x->data.F32[i] / map->binning->nXbin;
    7779        int yRuff = y->data.F32[i] / map->binning->nYbin;
    … …  
    7981        int bin = xRuff + yRuff*map->map->numCols;
    8082       
    81         psVector *vector = vectors->data[bin];
    82         vector->data.F32[vector->n] = f->data.F32[i];
     83        psVector *pixels = pixelSets->data[bin];
     84        pixels->data.S32[pixels->n] = i;
    8385        psVectorExtend (vector, 4, 1);
    8486    }
    8587
     88    // accumulate the x,y coords for each point to calculate the mean position.
    8689    int Nx = map->map->numCols;
    8790    int Ny = map->map->numRows;
    … …  
    8992        for (int ix = 0; ix < Nx; ix++) {
    9093   
    91             // select the vector
    92             psVector *vector = vectors->data[ix + iy*Nx];
     94            // select the vector.  use psImageUnbin command to convert xRuff,yRuff to
     95            // xFine,yFine
     96            psVector *pixels = pixelSets->data[ix + iy*Nx];
     97
     98            psVector xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     99            psVector yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     100            psVector fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     101
     102            for (int i = 0; i < pixels->n; i++) {
     103                bin = pixels->data.S32[i];
     104                xCell->data.F32[i] = x->data.F32[bin];
     105                yCell->data.F32[i] = y->data.F32[bin];
     106                fCell->data.F32[i] = f->data.F32[bin];
     107            }
    93108
    94109            // reset the stats to avoid contamination from the previous loop
    … …  
    96111
    97112            // get the value
    98             if (psVectorStats (map->stats, vector, NULL, NULL, 0)) {
     113            // XXX need to supply a mask and skip the masked pixels when calculating the centroid
     114            // this will not in general be properly weighted...
     115            if (psVectorStats (map->stats, fCell, NULL, NULL, 0)) {
    99116                mask->data.U8[iy][ix] = 0;
    100117                // XXX ensure only one option is selected, or save both position and width
    101118                map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options);
     119
     120                // calculate the mean position and save:
     121                psVectorStats (meanStat, xCell, NULL, NULL, 0);
     122                xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
     123                psVectorStats (meanStat, yCell, NULL, NULL, 0);
     124                yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
    102125            } else {
    103126                mask->data.U8[iy][ix] = 1;
    104127            }
     128           
     129            psFree (xCell);
     130            psFree (yCell);
     131            psFree (fCell);
    105132        }
    106133    }
     134
     135    // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
    107136
    108137    psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff);
    … …  
    112141    }
    113142
     143    // fit the valid pixels to (0,1,2) order polynomials, interpolate values to the pixel center
     144    // XXX I need to be careful about the pixel coordinates: center is 0,0 or 0.5, 0.5?
     145    psImagePixelInterpolateCenters (map->map, xCoord, yCoord, state, mask, 0xff);
     146    psFree (xCoord);
     147    psFree (yCoord);
     148
    114149    psImagePixelInterpolatePoor (map->map, state, mask, 0xff);
     150
    115151    return true;
    116152}

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

    r14725 r14726  
    1111 *  @author Eugene Magnier, IfA
    1212 *
    13  *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2007-09-02 20:29:50 $
     13 *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2007-09-03 20:29:07 $
    1515 *
    1616 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    163163                  break; }
    164164
     165                // XXX should I use 1 1D polynomial fitting all unmasked pixels in the 3x3 grid?
     166                // XXX that would automatically extend to regions where only 2 pixels are valid...
    165167              case PS_IMAGE_INTERPOLATE_LL: {
    166168                  // fit a plane to the 3 pixels at (0,1),(1,0),(1,1), extend to pixel at (0,0)
    … …  
    196198    return true;
    197199}
    198 
    199200   
     201// interpolate the good pixels to their true centers
     202bool psImagePixelInterpolateCenter (psImage *value, psImage *xCoord, psImage *yCoord, psImage *state, psImage *mask, psMaskType maskVal) {
     203
     204    assert (value->numCols == state->numCols);
     205    assert (value->numRows == state->numRows);
     206    assert (value->numCols == mask->numCols);
     207    assert (value->numRows == mask->numRows);
     208
     209    psImage *output = psImageAlloc (value->numCols, value->numRows, PS_TYPE_F32);
     210
     211    // allocate the vectors for the 2nd order fit below
     212    psVector *f  = psVectorAlloc (9, PS_TYPE_F32);
     213    // XXX if we add the weight above, include df
     214    // psVector *df = psVectorAlloc (9, PS_TYPE_F32);
     215    psVector *x  = psVectorAlloc (9, PS_TYPE_F32);
     216    psVector *y  = psVectorAlloc (9, PS_TYPE_F32);
     217
     218    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
     219    psPolynomial2D *poly2D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
     220    poly2D->mask[2][2] = 1;
     221    poly2D->mask[2][1] = 1;
     222    poly2D->mask[1][2] = 1;
     223
     224    // allocate a 2D polynomial to fit a plane to the valid neighbor pixels.
     225    psPolynomial2D *poly1D = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
     226    poly2D->mask[1][1] = 1;
     227
     228    for (int iy = 0; iy < state->numRows; iy++) {
     229        for (int ix = 0; ix < state->numCols; ix++) {
     230
     231            switch (state->data.S32[iy][ix]) {
     232              case PS_IMAGE_INTERPOLATE_GOOD2: {
     233                  // XXX is there a fit-image-region function?
     234                  int n = 0;
     235                  for (int jy = -1; jy <= +1; jy++) {
     236                      // skip invalid pixels
     237                      if (jy + iy < 0) { continue; }
     238                      if (jy + iy >= value->numRows) { continue; }
     239                      for (int jx = -1; jx <= +1; jx++) {
     240                          // skip invalid pixels
     241                          if (jx + ix < 0) { continue; }
     242                          if (jx + ix >= value->numCols) { continue; }
     243                          // skip masked pixels
     244                          if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
     245                          x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
     246                          y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
     247                          f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
     248                          // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
     249                          n++;
     250                      }
     251                  }
     252                  // set vector lengths here
     253                  x->n = n;
     254                  y->n = n;
     255                  f->n = n;
     256                  // df->n = n;
     257                  // psVectorFitPolynomial2D (poly, NULL, 0xff, f, df, x, y);
     258                  psVectorFitPolynomial2D (poly2D, NULL, 0xff, f, NULL, x, y);
     259                  // apply the fitted quadratic to get the poor pixel value
     260                  output->data.F32[iy][ix] = psPolynomial2DEval (poly2D, ix, iy);
     261                  break; }
     262
     263              case PS_IMAGE_INTERPOLATE_GOOD1: {
     264                  // XXX is there a fit-image-region function?
     265                  int n = 0;
     266                  for (int jy = -1; jy <= +1; jy++) {
     267                      // skip invalid pixels
     268                      if (jy + iy < 0) { continue; }
     269                      if (jy + iy >= value->numRows) { continue; }
     270                      for (int jx = -1; jx <= +1; jx++) {
     271                          // skip invalid pixels
     272                          if (jx + ix < 0) { continue; }
     273                          if (jx + ix >= value->numCols) { continue; }
     274                          // skip masked pixels
     275                          if (mask->data.PS_TYPE_MASK_DATA[iy+jy][ix+jx] & maskVal) { continue; }
     276                          x->data.F32[n] = xCoord->data.F32[iy+jy][ix+jx];
     277                          y->data.F32[n] = yCoord->data.F32[iy+jy][ix+jx];
     278                          f->data.F32[n] = value->data.F32[iy+jy][ix+jx];
     279                          // df->data.F32[n] = weight->data.F32[iy+jy][ix+jx];
     280                          n++;
     281                      }
     282                  }
     283                  // set vector lengths here
     284                  x->n = n;
     285                  y->n = n;
     286                  f->n = n;
     287                  // df->n = n;
     288                  // psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, df, x, y);
     289                  psVectorFitPolynomial2D (poly1D, NULL, 0xff, f, NULL, x, y);
     290                  // apply the fitted quadratic to get the poor pixel value
     291                  output->data.F32[iy][ix] = psPolynomial2DEval (poly1D, ix, iy);
     292                  break; }
     293
     294              case PS_IMAGE_INTERPOLATE_GOOD0: {
     295                  output->data.F32[iy][ix] = value->data.F32[iy+jy][ix+jx];
     296                  break; }
     297
     298              default:
     299                psAbort("impossible case in __func__");
     300            }
     301        }           
     302    }
     303
     304    for (int iy = 0; iy < value->numRows; iy++) {
     305        for (int ix = 0; ix < value->numCols; ix++) {
     306          if (mask->data.PS_TYPE_MASK_DATA[iy][ix] & maskVal) { continue; }
     307          value->data.F32[iy][ix] = output->data.F32[iy][ix];
     308        }
     309    }
     310
     311    psFree (x);
     312    psFree (y);
     313    psFree (f);
     314
     315    psFree (poly2D);
     316    psFree (poly1D);
     317
     318    psFree (output);
     319    return true;
     320}

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

    r14725 r14726  
    9292
    9393    // make a model for a well-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)
     94    {
     95        // function is defined over the range 0-1000, 0-1000
     96        psImageBinning *binning = psImageBinningAlloc();
     97        binning->nXfine = 1000;
     98        binning->nYfine = 1000;
     99        binning->nXruff = 5;
     100        binning->nYruff = 5;
     101
     102        // generate a grid of test data points
     103        psVector *x = psVectorAllocEmpty (100, PS_TYPE_F32);
     104        psVector *y = psVectorAllocEmpty (100, PS_TYPE_F32);
     105        psVector *f = psVectorAllocEmpty (100, PS_TYPE_F32);
     106
     107        for (int ix = 0; ix < 1000; ix += 10) {
     108            for (int iy = 0; iy < 1000; iy += 10) {
     109                x->data.F32[x->n] = ix;
     110                y->data.F32[y->n] = iy;
     111                f->data.F32[f->n] = C00 / (C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500) + 0.1);
     112                psVectorExtend (x, 100, 1);
     113                psVectorExtend (y, 100, 1);
     114                psVectorExtend (f, 100, 1);
     115            }
     116        }
     117
     118        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     119
     120        // scale defines both field and map image sizes (nXfine, nXruff)
     121        psImageMap *map = psImageMapAlloc (NULL, binning, stats);
     122
     123        // allocate a map, but we have no field image to supply
     124        // XXX this function needs to correct for the mean superpixel position
     125        // which is sampled...
     126        psImageMapGenerate (map, x, y, f, 0.1);
     127        psFree (binning);
     128
     129        fprintf (stderr, "nGood: %d\n", map->nGood);
     130        fprintf (stderr, "nPoor: %d\n", map->nPoor);
     131        fprintf (stderr, "nBad:  %d\n", map->nBad);
     132       
     133        SaveImage (NULL, map->map, "map.fits");
     134       
     135        psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32);
     136        for (int ix = 0; ix < 1000; ix++) {
     137            for (int iy = 0; iy < 1000; iy++) {
     138                field->data.F32[iy][ix] = C00 / (C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500) + 0.1);
     139            }
     140        }
     141        SaveImage (NULL, field, "field.fits");
     142
     143        // measure difference between model (map) and data (field)
     144        for (int ix = 0; ix < map->map->numCols; ix++) {
     145            for (int iy = 0; iy < map->map->numRows; iy++) {
     146             
     147                int xo = (ix + 0.5)*map->binning->nXbin;
     148                int yo = (iy + 0.5)*map->binning->nYbin;
     149                float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix];
     150                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]);
     151            }
     152        }
     153    }
    94154
    95155    // make a model for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)
Note: See TracChangeset for help on using the changeset viewer.