Changeset 14774

Timestamp:

Sep 7, 2007, 10:19:35 AM (19 years ago)

Author:

magnier

Message:

added additional functions to manage map construction

Location:

branches/eam_branch_20070830/psLib/src/imageops

Files:

• psImageMap.c (modified) (11 diffs)
• psImageMap.h (modified) (2 diffs)

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

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-03 20:29:07 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:19:35 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -20 +20 @@
 #include "psError.h"
 #include "psAbort.h"
+
+#include "psFits.h"
 #include "psAssert.h"
 #include "psRegion.h"
+#include "psFitsImage.h"
+
 #include "psVector.h"
 #include "psImage.h"
@@ -28 +32 @@
 #include "psImageMap.h"
 #include "psImagePixelInterpolate.h"
+#include "psImageUnbin.h"
 
 static void psImageMapFree(psImageMap *map) {
@@ -34 +39 @@
 
     psFree (map->map);
-    psFree (map->field);
+    // psFree (map->field); XXX ??? this should be freed here, but that causes an error...
     psFree (map->stats);
     psFree (map->binning);
@@ -55 +60 @@
 
     psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip (map->binning, map->field);
 
     return map;
 }
 
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff) {
+
+    assert (map);
+
+    map->binning->nXruff = nXruff;
+    map->binning->nYruff = nYruff;
+
+    psImageRecycle (map->map, nXruff, nYruff, PS_TYPE_F32);
+
+    psImageBinningSetScale (map->binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip (map->binning, map->field);
+
+    return true;
+}
+
 // generate a psImageMap (or NULL) with the given number of superpixels in X and Y
-bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac) {
+// this function returns an error if the output map has impossible holes
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
 
     psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK);
@@ -70 +92 @@
     // 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++) {
+    for (int i = 0; i < pixelSets->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;
+        int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
+        int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
 
         int bin = xRuff + yRuff*map->map->numCols;
+        assert (bin >= 0);
+        assert (bin < pixelSets->n);
         
         psVector *pixels = pixelSets->data[bin];
         pixels->data.S32[pixels->n] = i;
-        psVectorExtend (vector, 4, 1);
-    }
+        psVectorExtend (pixels, 4, 1);
+    }
+
+    // stats structure for getting the position centers
+    psStats *meanStat = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
 
     // accumulate the x,y coords for each point to calculate the mean position.
@@ -92 +118 @@
         for (int ix = 0; ix < Nx; ix++) {
     
-            // select the vector.  use psImageUnbin command to convert xRuff,yRuff to
-            // xFine,yFine
+            // pixel index for this cell
             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);
-
+            // storage vectors
+            psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+
+            // error vector, if needed
+            psVector *dfCell = NULL;
+            if (df) {
+                dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
+            }           
+
+            // collect data for this cell
             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];
+                int bin = pixels->data.S32[i];
+                // convert x,y in the fine image to the ruff image
+                xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
+                yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
+                fCell->data.F32[i]  = f->data.F32[bin];
+                if (df) {
+                    dfCell->data.F32[i] = df->data.F32[bin];
+                }
             }
 
@@ -113 +150 @@
             // 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)) {
+            if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
                 mask->data.U8[iy][ix] = 0;
                 // XXX ensure only one option is selected, or save both position and width
@@ -119 +156 @@
 
                 // calculate the mean position and save:
+                psStatsInit (meanStat);
                 psVectorStats (meanStat, xCell, NULL, NULL, 0);
                 xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
+                psStatsInit (meanStat);
                 psVectorStats (meanStat, yCell, NULL, NULL, 0);
                 yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options); 
@@ -130 +169 @@
             psFree (yCell);
             psFree (fCell);
-        }
-    }
-
-    // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
-
+            psFree (dfCell);
+        }
+    }
+    psFree (pixelSets);
+    psFree (meanStat);
+   // at this point, for each map pixel, we have (f,x,y), or the pixel is masked.
+
+    psFits *fits = NULL;
+
+    fits = psFitsOpen ("imageMap.raw.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    // did this analysis succeed?  (enough good or OK pixels?)
     psImage *state = psImagePixelInterpolateState (&map->nBad, &map->nPoor, mask, 0xff);
     map->nGood = mask->numCols * mask->numRows - map->nBad - map->nPoor;
     if (map->nBad > badFrac * mask->numCols * mask->numRows) {
+        psFree (xCoord);
+        psFree (yCoord);
+        psFree (mask);
         return false;
     }
@@ -143 +194 @@
     // 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);
+    psImagePixelInterpolateCenter (map->map, xCoord, yCoord, state, mask, 0xff);
     psFree (xCoord);
     psFree (yCoord);
 
+    fits = psFitsOpen ("imageMap.ref.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
     psImagePixelInterpolatePoor (map->map, state, mask, 0xff);
 
+    fits = psFitsOpen ("imageMap.fix.fits", "w");
+    psFitsWriteImage (fits, NULL, map->map, 0, NULL);
+    psFitsClose (fits);
+
+    psFree (state);
+    psFree (mask);
     return true;
 }
+
+// using the points given, generate a map with maximum resolution that yields only good and ok pixels
+bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
+    
+    int nXruff, nYruff;
+        
+    while (!psImageMapGenerate (map, x, y, f, df, badFrac)) {
+        // if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
+        // try again...  try to keep the aspect ratio.
+
+        // if both axes are at 1, give up
+        if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
+            return false;
+        }
+
+        // if one axis is at 1, decrement the other
+        if (map->binning->nXruff == 1) {
+            nXruff = map->binning->nXruff;
+            nYruff = map->binning->nYruff - 1;
+            psImageMapModifyScale (map, nXruff, nYruff);
+            continue;
+        }
+        if (map->binning->nYruff == 1) {
+            nYruff = map->binning->nYruff;
+            nXruff = map->binning->nXruff - 1;
+            psImageMapModifyScale (map, nXruff, nYruff);
+            continue;
+        }
+
+        // otherwise, decrement the larger axis, and set the smaller based
+        // on the aspect ratio
+        float aRatio = map->binning->nXruff / map->binning->nYruff;
+        if (map->binning->nXruff > map->binning->nYruff) {
+            nXruff = map->binning->nXruff - 1;
+            nYruff = (int)(0.5 + (nXruff / aRatio));
+        } else {
+            nYruff = map->binning->nYruff - 1;
+            nXruff = (int)(0.5 + (nYruff * aRatio));
+        }
+
+        psImageMapModifyScale (map, nXruff, nYruff);
+    }
+    return true;
+}
+
+double psImageMapEval (psImageMap *map, float x, float y) {
+
+    double result;
+
+    result = psImageUnbinPixel(x, y, map->map, map->binning);
+
+    return result; 
+}
+
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y) {
+
+    assert (x);
+    assert (y);
+    assert (x->n == y->n);
+    assert (map);
+
+    psVector *result = psVectorAlloc (x->n, PS_TYPE_F32);
+
+    for (int i = 0; i < x->n; i++) {
+        result->data.F32[i] = psImageUnbinPixel(x->data.F32[i], y->data.F32[i], map->map, map->binning);
+    }
+
+    return result; 
+}

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

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-02 02:03:58 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:19:35 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -35 +35 @@
 psImageMap *psImageMapAlloc(psImage *field, psImageBinning *binning, psStats *stats) PS_ATTR_MALLOC;
 
+bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff);
+
 // generate a psImageMap (or NULL) with the given number of superpixels in X and Y
-bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, float badFrac);
+bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
+
+// apply the psImageMap to the given coordinate (fine image pixels)
+double psImageMapEval (psImageMap *map, float x, float y);
+
+// apply the psImageMap to the given coordinate vectors (fine image pixels)
+psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y);
 
 /// @}