Changeset 10957

Timestamp:

Jan 7, 2007, 3:28:00 PM (20 years ago)

Author:

eugene

Message:

added local linear map

Location:

trunk/pswarp/src

Files:

• Makefile.am (modified) (1 diff)
• pswarp.h (modified) (2 diffs)
• pswarpMapGrid.c (modified) (2 diffs)
• pswarpMatchRange.c (added)
• pswarpTransformReadout.c (modified) (1 diff)
• pswarpTransformReadout_Opt.c (modified) (2 diffs)

trunk/pswarp/src/Makefile.am

@@ -18 +18 @@
 pswarpDefine.c             \
 pswarpErrorCodes.c         \
+pswarpMapGrid.c            \
+pswarpMatchRange.c         \
 pswarpParseCamera.c        \
 pswarpTransformReadout.c   \
+pswarpTransformReadout_Opt.c \
 pswarpVersion.c            
 

trunk/pswarp/src/pswarp.h

@@ -13 +13 @@
 # define PSWARP_RECIPE "PSWARP" // Name of the recipe to use
 
+// a single pswarpMap converts coordinates from one image to a second image 
+// the linear model is only valid over a limited range of pixels
 typedef struct {
-  double Xo, Xx, Xy;
-  double Yo, Yx, Yy;
+    double Xo, Xx, Xy;
+    double Yo, Yx, Yy;
+    int xo;
+    int yo;
 } pswarpMap;
 
+// the pswarpMapGrid carries a collection of pswarpMag structures representing the
+// local value of the pswarpMap at different locations in the image.
 typedef struct {
-  pswarpMap **grid;
-  int Nx, Ny;
-  int nXpix, nYpix;
+    pswarpMap ***maps;
+    int nXpts, nYpts;                   // number of x,y samples in the grid
+    int nXpix, nYpix;                   // x,y spacing in src image pixels of grid samples
+    int xMin,  yMin;                    // coordinate of first grid sample
 } pswarpMapGrid;
 
@@ -31 +38 @@
 void pswarpCleanup (pmConfig *config);
 bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);
+bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config);
 
+bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src);
+
+pswarpMap *pswarpMapAlloc ();
+pswarpMapGrid *pswarpMapGridAlloc (int Nx, int Ny);
+
+pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
+bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
+bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY);
+double pswarpMapGridMaxError (pswarpMapGrid *grid);
+bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY);
+bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy);

trunk/pswarp/src/pswarpMapGrid.c

@@ -1 +1 @@
 # include "pswarp.h"
 
-static pswarpMapFree (pswarpMap *map) {
+// construct a grid with superpixel spacing of nXpix, nYpix
+// XXX for the moment, ignore readout->cell->chip offsets
+pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) {
+
+    // split the difference of the remainder
+    int xMin = 0.5*(src->image->numCols % nXpix);
+    int yMin = 0.5*(src->image->numRows % nYpix);
+
+    pswarpMapGrid *grid = pswarpMapGridAlloc (src->image->numCols / nXpix, src->image->numRows / nYpix);
+
+    int ni = 0;
+    int nj = 0;
+    for (int i = xMin; i < src->image->numCols; i += nXpix, ni++) {
+        for (int j = yMin; j < src->image->numRows; j += nYpix, nj++) {
+            pswarpMapSetLocalModel (grid->maps[ni][nj], dest, src, i, j);
+        }
+    }
+            
+    grid->nXpix = nXpix;
+    grid->nYpix = nYpix;
+    grid->xMin = xMin;
+    grid->yMin = yMin;
+    return grid;
+}
+
+bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) {
+
+    *gridX = (ix - grid->xMin + 0.5*grid->nXpix) / grid->nXpix;
+    *gridY = (iy - grid->yMin + 0.5*grid->nYpix) / grid->nYpix;
+    return true;
+}
+
+bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY) {
+
+    *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix;
+    *nextY = gridY*grid->nYpix + grid->yMin + 0.5*grid->nYpix;
+    return true;
+}
+
+// measure the max error accumulated in appling one grid point to its neighbors
+double pswarpMapGridMaxError (pswarpMapGrid *grid) {
+
+    double xRaw, yRaw;
+    double xRef, yRef;
+    double maxError = 0;
+
+    for (int i = 0; i < grid->nXpts - 1; i++) {
+        for (int j = 0; j < grid->nYpts - 1; j++) {
+
+            // measure the output coordinates for the next grid position using the current grid map
+            // compare with the coordinates measured using the next grid map
+            pswarpMapApply (&xRaw, &yRaw, grid->maps[i][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo);
+            pswarpMapApply (&xRef, &yRef, grid->maps[i+1][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo);
+
+            double posError = hypot (xRaw-xRef, yRaw-yRef);
+            maxError = PS_MAX (maxError, posError);
+        }
+    }   
+    return maxError;
+}
+
+bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY) {
+
+    *outX = map->Xo + map->Xx*inX + map->Xy*inY;
+    *outY = map->Yo + map->Yx*inX + map->Yy*inY;
+
+    return true;
+}
+
+// determine the map for the given pixel from src to dest. pixel is in src coords
+bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy) {
+
+    pmCell *cell = NULL;
+
+    cell = src->parent;
+    pmChip *chipSrc = cell->parent;
+    pmFPA *fpaSrc = chipSrc->parent;
+
+    cell = dest->parent;
+    pmChip *chipDest = cell->parent;
+    pmFPA *fpaDest = chipDest->parent;
+
+    // XXX save these as static for speed?
+    psPlane *offset = psPlaneAlloc();
+
+    psPlane *FP = psPlaneAlloc();
+    psPlane *TP = psPlaneAlloc();
+    psSphere *sky = psSphereAlloc();
+
+    psPlane *V00 = psPlaneAlloc();
+    psPlane *V10 = psPlaneAlloc();
+    psPlane *V01 = psPlaneAlloc();
+
+    // XXX need to include readout->cell->chip offsets
+
+    // V(0,0) position
+    offset->x = ix;
+    offset->y = iy;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V00, chipDest->fromFPA, FP);
+    
+    // V(1,0) position
+    offset->x = ix + 1;
+    offset->y = iy;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V10, chipDest->fromFPA, FP);
+    
+    // V(0,1) position
+    offset->x = ix;
+    offset->y = iy + 1;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V01, chipDest->fromFPA, FP);
+
+    map->Xx = V10->x - V00->x;
+    map->Xy = V01->x - V00->x;
+    map->Xo = V00->x + map->Xx*ix + map->Xy*iy;
+
+    map->Yx = V10->y - V00->y;
+    map->Yy = V01->y - V00->y;
+    map->Yo = V00->y + map->Yx*ix + map->Yy*iy;
+ 
+    map->xo = ix;
+    map->yo = iy;
+
+    psFree (offset);
+    psFree (FP);
+    psFree (TP);
+    psFree (sky);
+
+    psFree (V00);
+    psFree (V10);
+    psFree (V01);
+
+    return true;
+}
+
+static void pswarpMapFree (pswarpMap *map) {
   return;
 }
@@ -13 +161 @@
 }
 
+static void pswarpMapGridFree (pswarpMapGrid *grid) {
+
+    if (grid == NULL) return;
+    if (grid->maps == NULL) return;
+
+    for (int i = 0; i < grid->nXpts; i++) {
+        for (int j = 0; j < grid->nYpts; j++) {
+            psFree (grid->maps[i][j]);
+        }
+        psFree (grid->maps[i]);
+    }
+    psFree (grid->maps);
+    return;
+}
+
+pswarpMapGrid *pswarpMapGridAlloc (int nXpts, int nYpts) {
+
+  pswarpMapGrid *grid = (pswarpMapGrid *) psAlloc (sizeof(pswarpMapGrid));
+  psMemSetDeallocator(grid, (psFreeFunc) pswarpMapGridFree);
+
+  grid->maps = psAlloc (nXpts*sizeof(void **));
+  for (int i = 0; i < nXpts; i++) {
+      grid->maps[i] = psAlloc (nYpts*sizeof(void *));
+      for (int j = 0; j < nYpts; j++) {
+          grid->maps[i][j] = pswarpMapAlloc();
+      }
+  }
+  grid->nXpts = nXpts;
+  grid->nYpts = nYpts;
+
+  grid->nXpix = 0;
+  grid->nYpix = 0;
   
+  return grid;
+}
+

trunk/pswarp/src/pswarpTransformReadout.c

@@ -35 +35 @@
         for (int x = 0; x < outNx; x++) {
             // Only transform those pixels requested
-            if (!region || (region && region->data.U8[y][x])) {
-                // Transform!
+            if (region && region->data.U8[y][x]) continue;
 
-                // XXX double check this 1/2 pixel offset
-                outPix->x = (double)x + 0.5;
-                outPix->y = (double)y + 0.5;
+            // XXX double check this 1/2 pixel offset
+            outPix->x = (double)x + 0.5;
+            outPix->y = (double)y + 0.5;
 
-                psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
-                psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
-                psDeproject (sky, TP, fpaOutput->toSky);
+            psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
+            psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
+            psDeproject (sky, TP, fpaOutput->toSky);
                 
-                psProject (TP, sky, fpaInput->toSky);
-                psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
-                psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
+            psProject (TP, sky, fpaInput->toSky);
+            psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
+            psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
 
-                // XXX get interpolation method from the recipe
-                outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
-                // modify zero and scale?
-
-                // outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR);
-                
-# if (0)
-                if (error) {
-                    // Error is actually the variance
-                    outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-                                                                                            detector->x,
-                                                                                            detector->y,
-                                                                                            mask, 1, NAN);
-                }
-                if (error) {
-                    outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
-                }
-# endif
-
-            } // Pixels of interest
-
+            // XXX get interpolation method from the recipe
+            outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
+            // outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR);
+            // modify zero and scale?
         }
-    } // Iterating over output pixels
-
+    }
     return true;
 }
 
-    
+# if (0)
+            if (error) {
+                // Error is actually the variance
+                outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+                                                                                        detector->x,
+                                                                                        detector->y,
+                                                                                        mask, 1, NAN);
+            }
+            if (error) {
+                outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
+            }
+# endif

trunk/pswarp/src/pswarpTransformReadout_Opt.c

@@ -1 +1 @@
 # include "pswarp.h"
 
-bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) {
+// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
+bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config) {
+
+    int minX, minY, maxX, maxY;
+    int gridX, gridY, nextGridX, nextGridY;
+    pswarpMap *map = NULL;
 
     // XXX this implementation currently ignores the use of the region 
     psImage *region = NULL;
-    pmCell *cell = NULL;
 
     // select the current recipe
-    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-
-    int outNx = output->image->numCols;
-    int outNy = output->image->numRows;
+    // psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
 
     psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
-    psPlane *outPix = psPlaneAlloc();   // Coordinates on the output detector
-
-    psPlane *FP = psPlaneAlloc();       // Coordinates on the focal plane
-    psPlane *TP = psPlaneAlloc();       // Coordinates on the tangent plane
-    psSphere *sky = psSphereAlloc();    // Coordinates on the sky
-
-    cell = input->parent;
-    pmChip *chipInput = cell->parent;
-    pmFPA *fpaInput = chipInput->parent;
-
-    cell = output->parent;
-    pmChip *chipOutput = cell->parent;
-    pmFPA *fpaOutput = chipOutput->parent;
-
+    psImage *inImage = input->image;
     psF32 **outData = output->image->data.F32;
-    psImage *inImage = input->image;
 
     // we might want to do the rectangular regions outside of the selection independently
@@ -37 +24 @@
     pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output);
 
+    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, 128, 128);
+
+    // XXX need to modify the grid based on this result and force the maxError < XXX
+    double maxError = pswarpMapGridMaxError (grid);
+    fprintf (stderr, "maximum error using this grid sampling: %f\n", maxError);
+
+    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
+    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+    map = grid->maps[gridX][gridY];
+
     // Iterate over the output image pixels
-    for (int y = 0; y < outNy; y++) {
-        for (int x = 0; x < outNx; x++) {
+    for (int y = minY; y < maxY; y++) {
+        if (y >= nextGridY) {
+            gridY ++;
+            pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+            map = grid->maps[gridX][gridY];
+        }
+        for (int x = minX; x < maxX; x++) {
+            if (x >= nextGridX) {
+                gridX ++;
+                pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+                map = grid->maps[gridX][gridY];
+            }
+
             // Only transform those pixels requested
-            if (!region || (region && region->data.U8[y][x])) {
-                // Transform!
+            if (region && region->data.U8[y][x]) continue;
 
-                // XXX double check this 1/2 pixel offset
-                outPix->x = (double)x + 0.5;
-                outPix->y = (double)y + 0.5;
+            // XXX double check this 1/2 pixel offset
+            // XXX subtract 0.5,0.5 from result?
+            pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
 
-                psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
-                psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
-                psDeproject (sky, TP, fpaOutput->toSky);
-                
-                psProject (TP, sky, fpaInput->toSky);
-                psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
-                psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
-
-                // XXX get interpolation method from the recipe
-                outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR);
-
-# if (0)
-                if (error) {
-                    // Error is actually the variance
-                    outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-                                                                                            detector->x,
-                                                                                            detector->y,
-                                                                                            mask, 1, NAN);
-                }
-# endif
-
-                outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
-# if (0)
-                if (error) {
-                    outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
-                }
-# endif
-
-            } // Pixels of interest
-
+            // XXX get interpolation method from the recipe
+            // XXX include mask
+            // XXX apply scale and offset?
+            outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
         }
-    } // Iterating over output pixels
-
+    }
     return true;
 }
 
+# if (0)
+    if (error) {
+        // Error is actually the variance
+        outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+                                                                                detector->x,
+                                                                                detector->y,
+                                                                                mask, 1, NAN);
+    }
+    if (error) {
+        outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
+    }
+# endif
+