Changeset 14726 for branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c

Timestamp:

Sep 3, 2007, 10:29:07 AM (19 years ago)

Author:

magnier

Message:

updates to image map code, tests

File:

• branches/eam_branch_20070830/psLib/src/imageops/psImageMap.c (modified) (6 diffs)

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

@@ -7 +7 @@
  *  @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 +63 @@
 
     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 +81 @@
         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 +92 @@
         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 +111 @@
 
             // 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 +141 @@
     }
 
+    // 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;
 }