Changeset 14768

Timestamp:

Sep 6, 2007, 5:27:02 PM (19 years ago)

Author:

Paul Price

Message:

Warp the source list as well as the pixels.

Location:

trunk/pswarp/src

Files:

• pswarpLoop.c (modified) (3 diffs)
• pswarpParseCamera.c (modified) (1 diff)
• pswarpTransformReadout_Opt.c (modified) (4 diffs)

trunk/pswarp/src/pswarpLoop.c

@@ -66 +66 @@
         pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT.MASK");
         pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT.WEIGHT");
+        pmFPAfileActivate(config->files, false, "PSWARP.OUTPUT.SOURCES");
         pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
         pmChip *chip;
@@ -89 +90 @@
         pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.MASK");
         pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.WEIGHT");
+        pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
     }
 
@@ -156 +158 @@
                 if (!readout->data_exists) {
                     continue;
+                }
+
+                // Copy the sources from the astrometry carrier to the input, so they can be accessed by
+                // pswarpTransformReadout
+                pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
+                psArray *sources = psMetadataLookupPtr(&mdok, astromRO->analysis,
+                                                       "PSPHOT.SOURCES"); // Sources from astrometry
+                if (sources) {
+                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY,
+                                     "Sources from input astrometry", sources);
                 }
 

trunk/pswarp/src/pswarpParseCamera.c

@@ -80 +80 @@
     }
 
+    if (astrom) {
+        pmFPAfile *outSources = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.SOURCES");
+        if (!outSources) {
+            psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.OUTPUT.SOURCES");
+            return false;
+        }
+        outSources->save = true;
+    }
+
+
     // Chip selection: turn on only the chips specified
     bool mdok;                          // Status of MD lookup

trunk/pswarp/src/pswarpTransformReadout_Opt.c

@@ -1 +1 @@
 # include "pswarp.h"
+
+#define SOURCE_ARRAY_BUFFER 100         // Size to grow the array of sources at a time
 
 // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
@@ -83 +85 @@
 
     // Iterate over the output image pixels (parent frame)
+    bool goodPixels = false;            // Any input pixels landing on the output image?
     for (int y = minY; y < maxY; y++) {
         if (y >= nextGridY) {
@@ -111 +114 @@
             if (inPix->y - inRow0 < 0) continue;
             if (inPix->y - inRow0 >= inImage->numRows) continue;
+
+            goodPixels = true;
 
             // XXX include mask
@@ -135 +140 @@
     }
 
-    output->data_exists = true;
-
     psFree(interp);
     psFree (inPix);
     psFree (grid);
+
+    // Transform sources
+    bool mdok;
+    psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source
+    if (goodPixels && mdok && inSources) {
+        pswarpMapGrid *sourceGrid = pswarpMapGridFromImage(output, input, nGridX, nGridY); // Grid for sources
+
+        psArray *outSources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, output->analysis,
+                                                                      "PSPHOT.SOURCES")); // Target sources
+        if (!outSources) {
+            outSources = psArrayAllocEmpty(SOURCE_ARRAY_BUFFER);
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY,
+                             "Warped sources", outSources);
+        }
+
+        for (int i = 0; i < inSources->n; i++) {
+            pmSource *source = inSources->data[i]; // Source of interest
+            pmModel *model = source->modelPSF; // Model for this source
+            float xIn, yIn;             // Coordinates of source
+            xIn = model->params->data.F32[PM_PAR_XPOS] - inCol0;
+            yIn = model->params->data.F32[PM_PAR_YPOS] - inRow0;
+            int xGrid, yGrid;           // Grid coordinates for local map
+            if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get grid coordinates for source at %f,%f\n",
+                        xIn, yIn);
+                psFree(outSources);
+                psFree(sourceGrid);
+                return false;
+            }
+            if (xGrid < 0 || xGrid >= sourceGrid->nXpts || yGrid < 0 || yGrid >= sourceGrid->nYpts) {
+                // It's not even on the grid
+                continue;
+            }
+
+            pswarpMap *map = sourceGrid->maps[xGrid][yGrid]; // Locally linear transformation
+            double xOut, yOut;          // Output coordinates
+            if (!pswarpMapApply(&xOut, &yOut, map, xIn, yIn)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to transform coordinates for source at %f,%f\n",
+                        xIn, yIn);
+                psFree(outSources);
+                psFree(sourceGrid);
+                return false;
+            }
+            xOut += outCol0;
+            yOut += outRow0;
+            if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) {
+                // It's not in the output image
+                continue;
+            }
+
+            // Generate the new source in the output frame
+            pmSource *new = pmSourceAlloc(); // New source
+            new->peak = pmPeakAlloc(xOut, yOut, source->peak->flux, PM_PEAK_LONE);
+            new->peak->flux = source->peak->flux;
+            new->type = PM_SOURCE_TYPE_STAR;
+            new->psfMag = source->psfMag;
+            new->errMag = source->errMag;
+            new->sky = source->sky;
+            new->skyErr = source->skyErr;
+            new->pixWeight = source->pixWeight;
+            new->modelPSF = pmModelAlloc(source->modelPSF->type);
+
+#if 0
+            // XXX Note that this will not set the correct axes
+            pmPSF_AxesToModel(new->modelPSF->params->data.F32,
+                              pmPSF_ModelToAxes(source->modelPSF->params->data.F32, 20.0));
+#endif
+
+            // Propagate the position erorrs
+            float dxIn, dyIn;           // Errors in input coordinates
+            dxIn = model->dparams->data.F32[PM_PAR_XPOS];
+            dyIn = model->dparams->data.F32[PM_PAR_YPOS];
+
+            float dxOut, dyOut;         // Errors in output coordinates
+            dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn));
+            dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn));
+
+            new->modelPSF->params->data.F32[PM_PAR_XPOS] = xOut;
+            new->modelPSF->params->data.F32[PM_PAR_YPOS] = yOut;
+            new->modelPSF->dparams->data.F32[PM_PAR_XPOS] = dxOut;
+            new->modelPSF->dparams->data.F32[PM_PAR_YPOS] = dyOut;
+
+            psArrayAdd(outSources, SOURCE_ARRAY_BUFFER, new);
+            psFree(new);                // Drop reference
+        }
+        psFree(sourceGrid);
+        psFree(outSources);             // Drop reference
+    }
+
+    output->data_exists = true;
+
     psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));