IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27437


Ignore:
Timestamp:
Mar 24, 2010, 1:59:47 PM (16 years ago)
Author:
Paul Price
Message:

Adding ability to get back to raw cell coordinates by providing a raw file (this presumably reads the image, which is unnecessary, but that's OK). Default input is now degrees, radians can be specified with the '-radians' command-line option. Default output is only defined coordinates which can be turned off with the '-all' command-line option.

Location:
branches/pap_stack/ppViz/src/ppCoord
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_stack/ppViz/src/ppCoord/ppCoord.h

    r26377 r27437  
    1010typedef struct {
    1111    psString astromName;                // Filename with astrometry
     12    psString rawName;                   // Filename with raw image (or NULL)
    1213    psString pixelsName;                // Filename with pixel coordinates
    1314    psString chipName;                  // Name of chip of interest
    1415    psString radecName;                 // Filename with sky coordinates
    1516    pmConfig *config;                   // Configuration
     17    bool radians;                       // RA,Dec are in radians?
     18    bool all;                           // Only all coordinates?
    1619} ppCoordData;
    1720

  • branches/pap_stack/ppViz/src/ppCoord/ppCoordArguments.c

    r26377 r27437  
    4747    psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
    4848    psMetadataAddStr(arguments, PS_LIST_TAIL, "-astrom", 0, "Filename with astrometry", NULL);
     49    psMetadataAddStr(arguments, PS_LIST_TAIL, "-raw", 0, "Filename with raw data", NULL);
    4950    psMetadataAddStr(arguments, PS_LIST_TAIL, "-pixels", 0, "Filename with pixel coordinates", NULL);
    5051    psMetadataAddStr(arguments, PS_LIST_TAIL, "-chip", 0, "Chip for pixel coordinates", NULL);
    51     psMetadataAddStr(arguments, PS_LIST_TAIL, "-radec", 0, "Filename with RA, Dec", NULL);
     52    psMetadataAddStr(arguments, PS_LIST_TAIL, "-radec", 0, "Filename with RA,Dec (default decimal degrees)", NULL);
     53    psMetadataAddBool(arguments, PS_LIST_TAIL, "-radians", 0, "RA,Dec in radians?", NULL);
     54    psMetadataAddBool(arguments, PS_LIST_TAIL, "-all", 0, "Output all coordinates?", NULL);
    5255    if (!psArgumentParse(arguments, &argc, argv) || argc != 1) {
    5356        usage(argv[0], arguments, data);
    5457    }
    5558
    56     bool mdok;                          // Status of MD lookup
    5759    data->astromName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-astrom"));
     60    data->rawName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-raw"));
    5861    data->pixelsName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-pixels"));
    5962    data->chipName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-chip"));
    60     data->radecName = psMemIncrRefCounter(psMetadataLookupStr(&mdok, arguments, "-radec"));
     63    data->radecName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-radec"));
     64    data->radians = psMetadataLookupBool(NULL, arguments, "-radians");
     65    data->all = psMetadataLookupBool(NULL, arguments, "-all");
    6166//    psMetadataAddStr(data->config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root name", data->outRoot);
    6267

  • branches/pap_stack/ppViz/src/ppCoord/ppCoordCamera.c

    r26377 r27437  
    3939    }
    4040
     41    if (data->rawName) {
     42        fileArguments("RAW", data->rawName, "Input raw image", data->config);
     43        pmFPAfile *raw = pmFPAfileDefineFromArgs(&status, data->config, "PPIMAGE.INPUT", "RAW"); // File
     44        if (!status || !raw) {
     45            psError(PS_ERR_IO, false, "Failed to build file from PPIMAGE.INPUT");
     46            return false;
     47        }
     48    }
     49
    4150#if 0
    4251    // Now the camera has been determined, we can read the recipe

  • branches/pap_stack/ppViz/src/ppCoord/ppCoordData.c

    r26377 r27437  
    1414{
    1515    psFree(data->astromName);
     16    psFree(data->rawName);
    1617    psFree(data->pixelsName);
    1718    psFree(data->chipName);
    … …  
    2829
    2930    data->astromName = NULL;
     31    data->rawName = NULL;
    3032    data->pixelsName = NULL;
    3133    data->chipName = NULL;

  • branches/pap_stack/ppViz/src/ppCoord/ppCoordLoop.c

    r26377 r27437  
    2020    }
    2121
     22    pmFPAfile *rawFile = data->rawName ? pmFPAfileSelectSingle(config->files, "PPIMAGE.INPUT", 0) :
     23        NULL; // File with raw image
     24
    2225    psArray *pixels = NULL, *radec = NULL; // Array of vectors with coordinates
    2326    psArray *radecOut = NULL;              // Output for sky coordinates
    … …  
    2528        pixels = psVectorsReadFromFile(data->pixelsName, "%f %f");
    2629        if (!pixels || pixels->n != 2) {
    27             psError(PS_ERR_UNKNOWN, false, "Unable to read pixel coordinates");
     30            psError(psErrorCodeLast(), false, "Unable to read pixel coordinates");
    2831            return false;
    2932        }
    … …  
    3235        radec = psVectorsReadFromFile(data->radecName, "%lf %lf");
    3336        if (!radec || radec->n != 2) {
    34             psError(PS_ERR_UNKNOWN, false, "Unable to read sky coordinates");
     37            psError(psErrorCodeLast(), false, "Unable to read sky coordinates");
    3538            return false;
    3639        }
    3740        psVector *ra = radec->data[0];  // RA coordinates
    3841        long num = ra->n;               // Number of coordinates
    39         radecOut = psArrayAlloc(3);
     42        radecOut = psArrayAlloc(4);
    4043        radecOut->data[0] = psArrayAlloc(num);
    4144        radecOut->data[1] = psVectorAlloc(num, PS_TYPE_F32);
    4245        radecOut->data[2] = psVectorAlloc(num, PS_TYPE_F32);
     46        radecOut->data[3] = psArrayAlloc(num);
    4347        psVectorInit(radecOut->data[1], NAN);
    4448        psVectorInit(radecOut->data[2], NAN);
    … …  
    6165    if (bilevelAstrometry) {
    6266        if (!pmAstromReadBilevelMosaic(astromFile->fpa, phu->header)) {
    63             psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
     67            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    6468            psFree(view);
    6569            return false;
    … …  
    7276            continue;
    7377        }
     78        if (rawFile) {
     79            pmChip *rawChip = pmFPAviewThisChip(view, rawFile->fpa);                   // Chip with raw data
     80            if (!rawChip || !rawChip->file_exists) {
     81                continue;
     82            }
     83        }
    7484        const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    7585        if (data->chipName && strcmp(chipName, data->chipName) != 0) {
    76                 continue;
     86            continue;
    7787        }
    7888        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    79             psError(PS_ERR_UNKNOWN, false, "Error loading data from files.");
     89            psError(psErrorCodeLast(), false, "Error loading data from files.");
    8090            return false;
    8191        }
    … …  
    8999        if (bilevelAstrometry) {
    90100            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    91                 psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel chip astrometry for input FPA.");
     101                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    92102                psFree(view);
    93103                return false;
    … …  
    97107            psWarning("Reading WCS astrometry for chip %s.", chipName);
    98108            if (!pmAstromReadWCS(astromFile->fpa, chip, hdu->header, 1.0)) {
    99                 psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
     109                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    100110                psFree(view);
    101111                return false;
    … …  
    120130                psDeproject(sky, tp, astromFile->fpa->toSky);
    121131
     132                if (!data->radians) {
     133                    sky->r *= M_PI / 180.0;
     134                    sky->d *= M_PI / 180.0;
     135                }
     136
    122137                fprintf(stdout, "%s %.3f %.3f --> %.10lf %.10lf\n", chipName, pix->x, pix->y, sky->r, sky->d);
    123138            }
    … …  
    129144
    130145        if (radec) {
     146            pmChip *rawChip = rawFile ? pmFPAviewThisChip(view, rawFile->fpa) : NULL; // Chip with raw
     147            psArray *cellBounds = NULL;                                               // Bounds of each cell
     148            psArray *cellNames = NULL;                                                // Names of each cell
     149            psVector *cellParityX = NULL, *cellParityY = NULL;                        // Parity for each cell
     150            psVector *cellX0 = NULL, *cellY0 = NULL;                                  // Offset for each cell
     151            psVector *cellBinX = NULL, *cellBinY = NULL;                              // Binning for each cell
     152            if (rawChip) {
     153                if (!rawChip->data_exists) {
     154                    // Not interested in this chip
     155                    continue;
     156                }
     157                int numCells = rawChip->cells->n; // Number of cells
     158
     159                cellBounds = psArrayAlloc(numCells);
     160                cellNames = psArrayAlloc(numCells);
     161                cellParityX = psVectorAlloc(numCells, PS_TYPE_S32);
     162                cellParityY = psVectorAlloc(numCells, PS_TYPE_S32);
     163                cellX0 = psVectorAlloc(numCells, PS_TYPE_S32);
     164                cellY0 = psVectorAlloc(numCells, PS_TYPE_S32);
     165                cellBinX = psVectorAlloc(numCells, PS_TYPE_S32);
     166                cellBinY = psVectorAlloc(numCells, PS_TYPE_S32);
     167
     168                pmCell *rawCell;        // Cell with raw data
     169                while ((rawCell = pmFPAviewNextCell(view, rawFile->fpa, 1))) {
     170                    if (!rawCell->data_exists) {
     171                        continue;
     172                    }
     173                    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     174                        psError(psErrorCodeLast(), false, "Error loading data from files.");
     175                        return false;
     176                    }
     177                    psRegion *region = pmCellExtent(rawCell); // Bounds of cell
     178                    if (!region) {
     179                        psError(psErrorCodeLast(), false, "Unable to determine extent of cell.");
     180                        return false;
     181                    }
     182                    cellBounds->data[view->cell] = region;
     183                    cellNames->data[view->cell] = psMemIncrRefCounter(psMetadataLookupStr(NULL,
     184                                                      rawCell->concepts, "CELL.NAME"));
     185                    cellParityX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XPARITY");
     186                    cellParityY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YPARITY");
     187                    cellX0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.X0");
     188                    cellY0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.Y0");
     189                    cellBinX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XBIN");
     190                    cellBinY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YBIN");
     191                    if (cellParityX->data.S32[view->cell] == 0 || cellParityY->data.S32[view->cell] == 0 ||
     192                        cellBinX->data.S32[view->cell] == 0 || cellBinY->data.S32[view->cell] == 0) {
     193                        psError(PM_ERR_CONCEPTS, true, "Concepts aren't set: %d %d %d %d %d %d\n",
     194                                cellParityX->data.S32[view->cell], cellParityY->data.S32[view->cell],
     195                                cellX0->data.S32[view->cell], cellY0->data.S32[view->cell],
     196                                cellBinX->data.S32[view->cell], cellBinY->data.S32[view->cell]);
     197                        return false;
     198                    }
     199
     200                    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     201                        psError(psErrorCodeLast(), false, "Error freeing data from files.");
     202                        return false;
     203                    }
     204                }
     205            }
     206
    131207            psVector *ra = radec->data[0], *dec = radec->data[1]; // Pixel coordinates
    132208            long num = ra->n;                                     // Number of coordinates
    … …  
    135211            int numRows = psMetadataLookupS32(NULL, hdu->header, "IMNAXIS2"); // Number of rows
    136212            if (numCols <= 0 || numRows <= 0) {
    137                 psError(PS_ERR_UNKNOWN, false, "Unable to read size of chip.");
     213                psError(psErrorCodeLast(), false, "Unable to read size of chip.");
    138214                return false;
    139215            }
    … …  
    147223            psVector *xPix = radecOut->data[1];   // x coordinate for pixels
    148224            psVector *yPix = radecOut->data[2];   // y coordinate for pixels
     225            psArray *cellPix = radecOut->data[3]; // Cell for pixels
    149226
    150227            for (long i = 0; i < num; i++) {
    151228                sky->r = ra->data.F64[i];
    152229                sky->d = dec->data.F64[i];
     230
     231                if (!data->radians) {
     232                    sky->r *= M_PI / 180.0;
     233                    sky->d *= M_PI / 180.0;
     234                }
    153235
    154236                psProject(tp, sky, astromFile->fpa->toSky);
    … …  
    162244                }
    163245
     246                if (rawChip) {
     247                    psString cellName = NULL;         // Name of cell
     248                    int numCells = rawChip->cells->n; // Number of cells
     249                    for (int i = 0; i < numCells && !cellName; i++) {
     250                        psRegion *region = cellBounds->data[i]; // Bounds of cell
     251                        if (x >= region->x0 && x < region->x1 && y >= region->y0 && y < region->y1) {
     252                            cellName = psStringCopy(cellNames->data[i]);
     253                            // Transform chip coordinates to cell coordinates
     254                            x = (x - cellX0->data.S32[i]) /
     255                                (float)(cellParityX->data.S32[i] * cellBinX->data.S32[i]);
     256                            y = (y - cellY0->data.S32[i]) /
     257                                (float)(cellParityY->data.S32[i] * cellBinY->data.S32[i]);
     258                        }
     259                    }
     260                    cellPix->data[i] = cellName;
     261                }
     262
    164263                chipPix->data[i] = psStringCopy(chipName);
    165264                xPix->data.F32[i] = x;
    … …  
    170269            psFree(tp);
    171270            psFree(sky);
     271            psFree(cellNames);
     272            psFree(cellBounds);
     273            psFree(cellParityX);
     274            psFree(cellParityY);
     275            psFree(cellBinX);
     276            psFree(cellBinY);
     277            psFree(cellX0);
     278            psFree(cellY0);
    172279        }
    173280
    174281        // Chip
    175282        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    176             psError(PS_ERR_UNKNOWN, false, "Error saving data to files.");
     283            psError(psErrorCodeLast(), false, "Error saving data to files.");
    177284            return false;
    178285        }
    … …  
    180287    // FPA
    181288    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    182         psError(PS_ERR_UNKNOWN, false, "Error saving data to files.");
     289        psError(psErrorCodeLast(), false, "Error saving data to files.");
    183290        return false;
    184291    }
    … …  
    188295        psVector *xPix = radecOut->data[1];   // x coordinate for pixels
    189296        psVector *yPix = radecOut->data[2];   // y coordinate for pixels
     297        psArray *cellPix = radecOut->data[3]; // Cell for pixels
    190298        psVector *ra = radec->data[0];        // RA coordinate
    191299        psVector *dec = radec->data[1];       // Dec coordinate
    … …  
    193301        for (long i = 0; i < chipPix->n; i++) {
    194302            const char *chipName = chipPix->data[i]; // Name of chip
    195             fprintf(stdout, "%.10lf %.10lf --> %.3f %.3f %s\n",
     303            const char *cellName = cellPix->data[i]; // Name of cell, or NULL
     304            if (!data->all && (!isfinite(xPix->data.F32[i]) || !isfinite(yPix->data.F32[i]) ||
     305                               !chipName || (rawFile && !cellName))) {
     306                continue;
     307            }
     308            fprintf(stdout, "%.10lf %.10lf --> %.3f %.3f %s%s%s\n",
    196309                    ra->data.F64[i], dec->data.F64[i], xPix->data.F32[i], yPix->data.F32[i],
    197                     chipName ? chipName : "UNKNOWN");
     310                    chipName ? chipName : "UNKNOWN",
     311                    rawFile ? " " : "",
     312                    rawFile && cellName ? cellName : (rawFile ? "UNKNOWN" : ""));
    198313        }
    199314    }
Note: See TracChangeset for help on using the changeset viewer.