Changeset 27978

Timestamp:

May 17, 2010, 12:20:57 PM (16 years ago)

Author:

Paul Price

Message:

Adding support for Magic clusters list: file stars with un-commented 'RA DEC', so need to parse it a bit differently. Also hard-coded the units of the streaks and clusters (radians for streaks, degrees for clusters).

Location:

trunk/ppViz/src/ppCoord

Files:

• ppCoord.h (modified) (1 diff)
• ppCoordArguments.c (modified) (2 diffs)
• ppCoordLoop.c (modified) (8 diffs)

trunk/ppViz/src/ppCoord/ppCoord.h

@@ -17 +17 @@
     psString radecName;                 // Filename with sky coordinates
     psString streaksName;               // Filename with streaks (sky coordinates)
+    psString clustersName;              // Filename with clusters (sky coordinates)
     pmConfig *config;                   // Configuration
     bool radians;                       // RA,Dec are in radians?

trunk/ppViz/src/ppCoord/ppCoordArguments.c

@@ -52 +52 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-radec", 0, "Filename with RA,Dec (default decimal degrees)", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-streaks", 0, "Filename with streaks", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-clusters", 0, "Filename with clusters", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-radians", 0, "RA,Dec in radians?", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-all", 0, "Output all coordinates?", NULL);
@@ -67 +68 @@
     data->radecName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-radec"));
     data->streaksName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-streaks"));
+    data->clustersName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-clusters"));
     data->radians = psMetadataLookupBool(NULL, arguments, "-radians");
     data->all = psMetadataLookupBool(NULL, arguments, "-all");

trunk/ppViz/src/ppCoord/ppCoordLoop.c

@@ -80 +80 @@
         NULL; // File with raw image
 
-    psArray *pixels = NULL, *radec = NULL, *streaks = NULL; // Array of vectors with coordinates
-    psArray *radecOut = NULL, *streaksOut = NULL;           // Output for sky coordinates
+    psArray *pixels = NULL, *radec = NULL, *streaks = NULL, *clusters = NULL; // Array of coordinate vectors
+    psArray *radecOut = NULL, *streaksOut = NULL, *clustersOut = NULL;        // Output for sky coordinates
     if (data->pixelsName) {
         pixels = psVectorsReadFromFile(data->pixelsName, "%f %f");
@@ -145 +145 @@
         psVectorInit(streaksOut->data[6], NAN);
         psVectorInit(streaksOut->data[7], NAN);
+    }
+
+    if (data->clustersName) {
+        psString file = psSlurpFilename(data->clustersName); // Contents of clusters file
+        if (!file) {
+            psError(psErrorCodeLast(), false, "Unable to read clusters file %s", data->clustersName);
+            return false;
+        }
+        psArray *lines = psStringSplitArray(file, "\n", false); // Lines of clusters
+        psFree(file);
+
+        int num = lines->n - 1;         // Number of clusters
+        clusters = psArrayAlloc(2);
+        psVector *ra = clusters->data[0] = psVectorAlloc(num, PS_TYPE_F64);
+        psVector *dec = clusters->data[1] = psVectorAlloc(num, PS_TYPE_F64);
+
+        // Skip the first line
+        for (int i = 1; i < lines->n; i++) {
+            const char *line = lines->data[i]; // Line of interest
+            if (!sscanf(line, "%lf %lf", &ra->data.F64[i-1], &dec->data.F64[i-1]) != 2) {
+                psError(PS_ERR_IO, true, "Unable to read line %d of %s", i, data->clustersName);
+                return false;
+            }
+        }
+        psFree(lines);
+
+        clustersOut = psArrayAlloc(4);
+        clustersOut->data[0] = psArrayAlloc(num);
+        clustersOut->data[1] = psVectorAlloc(num, PS_TYPE_F32);
+        clustersOut->data[2] = psVectorAlloc(num, PS_TYPE_F32);
+        clustersOut->data[3] = psArrayAlloc(num);
+        psVectorInit(clustersOut->data[1], NAN);
+        psVectorInit(clustersOut->data[2], NAN);
     }
 
@@ -369 +402 @@
                 float x1, y1;   // Coordinates of point 1
                 coordSky2Chip(&x1, &y1, ra1->data.F64[i], dec1->data.F64[i],
-                              data->radians, astromFile->fpa, chip);
+                              true, astromFile->fpa, chip);
                 if ((x1 >= 0 && x1 < numCols && y1 >= 0 && y1 < numRows)) {
                     if (rawChip) {
@@ -384 +417 @@
                 float x2, y2;   // Coordinates of point 2
                 coordSky2Chip(&x2, &y2, ra2->data.F64[i], dec2->data.F64[i],
-                              data->radians, astromFile->fpa, chip);
+                              true, astromFile->fpa, chip);
                 if ((x2 >= 0 && x2 < numCols && y2 >= 0 && y2 < numRows)) {
                     if (rawChip) {
@@ -396 +429 @@
                     yPix2->data.F32[i] = y2;
                 }
+            }
+        }
+
+        if (clusters) {
+            psVector *ra = clusters->data[0], *dec = clusters->data[1]; // Pixel coordinates
+            long num = ra->n;                                     // Number of coordinates
+
+            int numCols = psMetadataLookupS32(NULL, hdu->header, "IMNAXIS1"); // Number of columns
+            int numRows = psMetadataLookupS32(NULL, hdu->header, "IMNAXIS2"); // Number of rows
+            if (numCols <= 0 || numRows <= 0) {
+                psError(psErrorCodeLast(), false, "Unable to read size of chip.");
+                return false;
+            }
+
+            psArray *chipPix = clustersOut->data[0]; // Chip for pixels
+            psVector *xPix = clustersOut->data[1];   // x coordinate for pixels
+            psVector *yPix = clustersOut->data[2];   // y coordinate for pixels
+            psArray *cellPix = clustersOut->data[3]; // Cell for pixels
+
+            for (long i = 0; i < num; i++) {
+                float x, y;             // Pixel coordinates
+                coordSky2Chip(&x, &y, ra->data.F64[i], dec->data.F64[i],
+                                false, astromFile->fpa, chip);
+                if ((x < 0 || x > numCols || y < 0 || y > numRows)) {
+                    // Not on this chip
+                    continue;
+                }
+
+                if (rawChip) {
+                    psString cellName = NULL; // Name of cell
+                    coordChip2Cell(&cellName, &x, &y, x, y, cellNames, cellBounds,
+                                   cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY);
+                    cellPix->data[i] = cellName;
+                }
+
+                chipPix->data[i] = psStringCopy(chipName);
+                xPix->data.F32[i] = x;
+                yPix->data.F32[i] = y;
             }
         }
@@ -482 +553 @@
                         data->ds9color);
             } else {
-                fprintf(stdout, "%.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %s%s%s %.3f %.3f %s%s%s\n",
+                fprintf(stdout, "Streak %.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %s%s%s %.3f %.3f %s%s%s\n",
                         ra1->data.F64[i], dec1->data.F64[i],
                         ra2->data.F64[i], dec2->data.F64[i],
@@ -498 +569 @@
     }
 
+    if (clustersOut) {
+        psArray *chipPix = clustersOut->data[0]; // Chip for pixels
+        psVector *xPix = clustersOut->data[1];   // x coordinate for pixels
+        psVector *yPix = clustersOut->data[2];   // y coordinate for pixels
+        psArray *cellPix = clustersOut->data[3]; // Cell for pixels
+        psVector *ra = clusters->data[0];        // RA coordinate
+        psVector *dec = clusters->data[1];       // Dec coordinate
+
+        for (long i = 0; i < chipPix->n; i++) {
+            const char *chipName = chipPix->data[i]; // Name of chip
+            const char *cellName = cellPix->data[i]; // Name of cell, or NULL
+            if (!data->all && (!isfinite(xPix->data.F32[i]) || !isfinite(yPix->data.F32[i]) ||
+                               !chipName || (rawFile && !cellName))) {
+                continue;
+            }
+            if (!rawFile && data->ds9) {
+                // Region file is only appropriate if we're not mapping all the way back to cell coordinates
+                fprintf(data->ds9, "image;circle(%f,%f,%f) # color=%s\n",
+                        xPix->data.F32[i], yPix->data.F32[i], data->ds9radius, data->ds9color);
+            } else {
+                fprintf(stdout, "Cluster %.10lf %.10lf --> %.3f %.3f %s%s%s\n",
+                        ra->data.F64[i], dec->data.F64[i], xPix->data.F32[i], yPix->data.F32[i],
+                        chipName ? chipName : "UNKNOWN",
+                        rawFile ? " " : "",
+                        rawFile && cellName ? cellName : (rawFile ? "UNKNOWN" : ""));
+            }
+        }
+    }
+
     psFree(pixels);
     psFree(radec);
@@ -503 +603 @@
     psFree(streaks);
     psFree(streaksOut);
+    psFree(clusters);
+    psFree(clustersOut);
 
     return true;