Changeset 10957 for trunk/pswarp/src/pswarpTransformReadout_Opt.c

Timestamp:

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

Author:

eugene

Message:

added local linear map

File:

• trunk/pswarp/src/pswarpTransformReadout_Opt.c (modified) (2 diffs)

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
+