IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 5, 2010, 6:26:48 PM (16 years ago)
Author:
Paul Price
Message:

Reorganising coordinate conversions to make it easy to calculate raw coordinates for streaks.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppViz/src/ppCoord/ppCoordLoop.c

    r27609 r27613  
    99#include "ppCoord.h"
    1010
     11// Convert sky coordinates to chip coordinates
     12static void coordSky2Chip(float *x, float *y, // Chip coordinates (output)
     13                          double ra, double dec, // Sky coordinates
     14                          bool radians,          // Coordinates are in radians?
     15                          const pmFPA *astromFPA,  // Astrometry FPA
     16                          const pmChip *astromChip // Astrometry chip
     17    )
     18{
     19    psPlane pix;                        // Pixel coordinates on chip
     20    psPlane fp;                         // Focal plane coordinates
     21    psPlane tp;                         // Tangent plane coordinates
     22    psSphere sky;                       // Sky coordinates
     23
     24    sky.r = ra;
     25    sky.d = dec;
     26
     27    if (!radians) {
     28        sky.r *= M_PI / 180.0;
     29        sky.d *= M_PI / 180.0;
     30    }
     31
     32    psProject(&tp, &sky, astromFPA->toSky);
     33    psPlaneTransformApply(&fp, astromFPA->fromTPA, &tp);
     34    psPlaneTransformApply(&pix, astromChip->fromFPA, &fp);
     35
     36    *x = pix.x;
     37    *y = pix.y;
     38
     39    return;
     40}
     41
     42// Convert chip coordinates to cell coordinates
     43static void coordChip2Cell(psString *cellName,       // Cell name (output)
     44                           float *xCell, float *yCell, // Pixel coordinates (output)
     45                           float xChip, float yChip, // Sky coordinates
     46                           const psArray *cellNames, // Names of cells
     47                           const psArray *cellBounds,                                // Bounds of cells
     48                           const psVector *cellX0, const psVector *cellY0, // Cell offsets
     49                           const psVector *cellParityX, const psVector *cellParityY, // Cell parities
     50                           const psVector *cellBinX, const psVector *cellBinY        // Cell binnings
     51
     52
     53    )
     54{
     55    int numCells = cellNames->n; // Number of cells
     56    for (int i = 0; i < numCells && !cellName; i++) {
     57        psRegion *region = cellBounds->data[i]; // Bounds of cell
     58        if (xChip >= region->x0 && xChip < region->x1 && yChip >= region->y0 && yChip < region->y1) {
     59            *cellName = psStringCopy(cellNames->data[i]);
     60            // Transform chip coordinates to cell coordinates
     61            *xCell = (xChip - cellX0->data.S32[i]) / (float)(cellParityX->data.S32[i] * cellBinX->data.S32[i]);
     62            *yCell = (yChip - cellY0->data.S32[i]) / (float)(cellParityY->data.S32[i] * cellBinY->data.S32[i]);
     63        }
     64    }
     65
     66    return;
     67}
    1168
    1269bool ppCoordLoop(ppCoordData *data // Run-time data
    … …  
    75132            }
    76133        }
    77         streaksOut = psArrayAlloc(4);
    78         streaksOut->data[0] = psVectorAlloc(numStreaks, PS_TYPE_F32);
    79         streaksOut->data[1] = psVectorAlloc(numStreaks, PS_TYPE_F32);
     134        streaksOut = psArrayAlloc(8);
     135        streaksOut->data[0] = psArrayAlloc(numStreaks);
     136        streaksOut->data[1] = psArrayAlloc(numStreaks);
    80137        streaksOut->data[2] = psVectorAlloc(numStreaks, PS_TYPE_F32);
    81138        streaksOut->data[3] = psVectorAlloc(numStreaks, PS_TYPE_F32);
    82         psVectorInit(streaksOut->data[0], NAN);
    83         psVectorInit(streaksOut->data[1], NAN);
    84         psVectorInit(streaksOut->data[2], NAN);
     139        streaksOut->data[4] = psArrayAlloc(numStreaks);
     140        streaksOut->data[5] = psArrayAlloc(numStreaks);
     141        streaksOut->data[6] = psVectorAlloc(numStreaks, PS_TYPE_F32);
     142        streaksOut->data[7] = psVectorAlloc(numStreaks, PS_TYPE_F32);
    85143        psVectorInit(streaksOut->data[3], NAN);
     144        psVectorInit(streaksOut->data[4], NAN);
     145        psVectorInit(streaksOut->data[6], NAN);
     146        psVectorInit(streaksOut->data[7], NAN);
    86147    }
    87148
    … …  
    183244        }
    184245
     246
     247        pmChip *rawChip = rawFile ? pmFPAviewThisChip(view, rawFile->fpa) : NULL; // Chip with raw
     248        psArray *cellBounds = NULL;                                               // Bounds of each cell
     249        psArray *cellNames = NULL;                                                // Names of each cell
     250        psVector *cellParityX = NULL, *cellParityY = NULL;                        // Parity for each cell
     251        psVector *cellX0 = NULL, *cellY0 = NULL;                                  // Offset for each cell
     252        psVector *cellBinX = NULL, *cellBinY = NULL;                              // Binning for each cell
     253        if (rawChip) {
     254            if (!rawChip->data_exists) {
     255                // Not interested in this chip
     256                continue;
     257            }
     258            int numCells = rawChip->cells->n; // Number of cells
     259
     260            cellBounds = psArrayAlloc(numCells);
     261            cellNames = psArrayAlloc(numCells);
     262            cellParityX = psVectorAlloc(numCells, PS_TYPE_S32);
     263            cellParityY = psVectorAlloc(numCells, PS_TYPE_S32);
     264            cellX0 = psVectorAlloc(numCells, PS_TYPE_S32);
     265            cellY0 = psVectorAlloc(numCells, PS_TYPE_S32);
     266            cellBinX = psVectorAlloc(numCells, PS_TYPE_S32);
     267            cellBinY = psVectorAlloc(numCells, PS_TYPE_S32);
     268
     269            pmCell *rawCell;        // Cell with raw data
     270            while ((rawCell = pmFPAviewNextCell(view, rawFile->fpa, 1))) {
     271                if (!rawCell->data_exists) {
     272                    continue;
     273                }
     274                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     275                    psError(psErrorCodeLast(), false, "Error loading data from files.");
     276                    return false;
     277                }
     278                psRegion *region = pmCellExtent(rawCell); // Bounds of cell
     279                if (!region) {
     280                    psError(psErrorCodeLast(), false, "Unable to determine extent of cell.");
     281                    return false;
     282                }
     283                cellBounds->data[view->cell] = region;
     284                cellNames->data[view->cell] = psMemIncrRefCounter(psMetadataLookupStr(NULL, rawCell->concepts, "CELL.NAME"));
     285                cellParityX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XPARITY");
     286                cellParityY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YPARITY");
     287                cellX0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.X0");
     288                cellY0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.Y0");
     289                cellBinX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XBIN");
     290                cellBinY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YBIN");
     291                if (cellParityX->data.S32[view->cell] == 0 || cellParityY->data.S32[view->cell] == 0 ||
     292                    cellBinX->data.S32[view->cell] == 0 || cellBinY->data.S32[view->cell] == 0) {
     293                    psError(PM_ERR_CONCEPTS, true, "Concepts aren't set: %d %d %d %d %d %d\n",
     294                            cellParityX->data.S32[view->cell], cellParityY->data.S32[view->cell],
     295                            cellX0->data.S32[view->cell], cellY0->data.S32[view->cell],
     296                            cellBinX->data.S32[view->cell], cellBinY->data.S32[view->cell]);
     297                    return false;
     298                }
     299
     300                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     301                    psError(psErrorCodeLast(), false, "Error freeing data from files.");
     302                    return false;
     303                }
     304            }
     305        }
     306
    185307        if (radec) {
    186             pmChip *rawChip = rawFile ? pmFPAviewThisChip(view, rawFile->fpa) : NULL; // Chip with raw
    187             psArray *cellBounds = NULL;                                               // Bounds of each cell
    188             psArray *cellNames = NULL;                                                // Names of each cell
    189             psVector *cellParityX = NULL, *cellParityY = NULL;                        // Parity for each cell
    190             psVector *cellX0 = NULL, *cellY0 = NULL;                                  // Offset for each cell
    191             psVector *cellBinX = NULL, *cellBinY = NULL;                              // Binning for each cell
    192             if (rawChip) {
    193                 if (!rawChip->data_exists) {
    194                     // Not interested in this chip
    195                     continue;
    196                 }
    197                 int numCells = rawChip->cells->n; // Number of cells
    198 
    199                 cellBounds = psArrayAlloc(numCells);
    200                 cellNames = psArrayAlloc(numCells);
    201                 cellParityX = psVectorAlloc(numCells, PS_TYPE_S32);
    202                 cellParityY = psVectorAlloc(numCells, PS_TYPE_S32);
    203                 cellX0 = psVectorAlloc(numCells, PS_TYPE_S32);
    204                 cellY0 = psVectorAlloc(numCells, PS_TYPE_S32);
    205                 cellBinX = psVectorAlloc(numCells, PS_TYPE_S32);
    206                 cellBinY = psVectorAlloc(numCells, PS_TYPE_S32);
    207 
    208                 pmCell *rawCell;        // Cell with raw data
    209                 while ((rawCell = pmFPAviewNextCell(view, rawFile->fpa, 1))) {
    210                     if (!rawCell->data_exists) {
    211                         continue;
    212                     }
    213                     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    214                         psError(psErrorCodeLast(), false, "Error loading data from files.");
    215                         return false;
    216                     }
    217                     psRegion *region = pmCellExtent(rawCell); // Bounds of cell
    218                     if (!region) {
    219                         psError(psErrorCodeLast(), false, "Unable to determine extent of cell.");
    220                         return false;
    221                     }
    222                     cellBounds->data[view->cell] = region;
    223                     cellNames->data[view->cell] = psMemIncrRefCounter(psMetadataLookupStr(NULL,
    224                                                       rawCell->concepts, "CELL.NAME"));
    225                     cellParityX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XPARITY");
    226                     cellParityY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YPARITY");
    227                     cellX0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.X0");
    228                     cellY0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.Y0");
    229                     cellBinX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XBIN");
    230                     cellBinY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YBIN");
    231                     if (cellParityX->data.S32[view->cell] == 0 || cellParityY->data.S32[view->cell] == 0 ||
    232                         cellBinX->data.S32[view->cell] == 0 || cellBinY->data.S32[view->cell] == 0) {
    233                         psError(PM_ERR_CONCEPTS, true, "Concepts aren't set: %d %d %d %d %d %d\n",
    234                                 cellParityX->data.S32[view->cell], cellParityY->data.S32[view->cell],
    235                                 cellX0->data.S32[view->cell], cellY0->data.S32[view->cell],
    236                                 cellBinX->data.S32[view->cell], cellBinY->data.S32[view->cell]);
    237                         return false;
    238                     }
    239 
    240                     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    241                         psError(psErrorCodeLast(), false, "Error freeing data from files.");
    242                         return false;
    243                     }
    244                 }
    245             }
    246 
    247308            psVector *ra = radec->data[0], *dec = radec->data[1]; // Pixel coordinates
    248309            long num = ra->n;                                     // Number of coordinates
    … …  
    255316            }
    256317
    257             psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
    258             psPlane *fp = psPlaneAlloc();    // Focal plane coordinates
    259             psPlane *tp = psPlaneAlloc();    // Tangent plane coordinates
    260             psSphere *sky = psSphereAlloc(); // Sky coordinates
    261 
    262318            psArray *chipPix = radecOut->data[0]; // Chip for pixels
    263319            psVector *xPix = radecOut->data[1];   // x coordinate for pixels
    … …  
    266322
    267323            for (long i = 0; i < num; i++) {
    268                 sky->r = ra->data.F64[i];
    269                 sky->d = dec->data.F64[i];
    270 
    271                 if (!data->radians) {
    272                     sky->r *= M_PI / 180.0;
    273                     sky->d *= M_PI / 180.0;
    274                 }
    275 
    276                 psProject(tp, sky, astromFile->fpa->toSky);
    277                 psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp);
    278                 psPlaneTransformApply(pix, chip->fromFPA, fp);
    279 
    280                 float x = pix->x, y = pix->y; // Pixel coordinates
     324                float x, y;             // Pixel coordinates
     325                coordSky2Chip(&x, &y, ra->data.F64[i], dec->data.F64[i],
     326                                data->radians, astromFile->fpa, chip);
    281327                if ((x < 0 || x > numCols || y < 0 || y > numRows)) {
    282328                    // Not on this chip
    … …  
    285331
    286332                if (rawChip) {
    287                     psString cellName = NULL;         // Name of cell
    288                     int numCells = rawChip->cells->n; // Number of cells
    289                     for (int i = 0; i < numCells && !cellName; i++) {
    290                         psRegion *region = cellBounds->data[i]; // Bounds of cell
    291                         if (x >= region->x0 && x < region->x1 && y >= region->y0 && y < region->y1) {
    292                             cellName = psStringCopy(cellNames->data[i]);
    293                             // Transform chip coordinates to cell coordinates
    294                             x = (x - cellX0->data.S32[i]) /
    295                                 (float)(cellParityX->data.S32[i] * cellBinX->data.S32[i]);
    296                             y = (y - cellY0->data.S32[i]) /
    297                                 (float)(cellParityY->data.S32[i] * cellBinY->data.S32[i]);
    298                         }
    299                     }
     333                    psString cellName = NULL; // Name of cell
     334                    coordChip2Cell(&cellName, &x, &y, x, y, cellNames, cellBounds,
     335                                   cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY);
    300336                    cellPix->data[i] = cellName;
    301337                }
    … …  
    305341                yPix->data.F32[i] = y;
    306342            }
    307             psFree(pix);
    308             psFree(fp);
    309             psFree(tp);
    310             psFree(sky);
    311             psFree(cellNames);
    312             psFree(cellBounds);
    313             psFree(cellParityX);
    314             psFree(cellParityY);
    315             psFree(cellBinX);
    316             psFree(cellBinY);
    317             psFree(cellX0);
    318             psFree(cellY0);
    319         }
    320 
    321         if (streaks && data->chipName && !rawFile) {
    322             psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1
    323             psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1
    324             psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2
    325             psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2
     343        }
     344
     345        if (streaks) {
     346            psArray *chipPix1 = streaksOut->data[0]; // Chip for point 1
     347            psArray *cellPix1 = streaksOut->data[1]; // Cell for point 1
     348            psVector *xPix1 = streaksOut->data[2]; // x coordinate for point 1
     349            psVector *yPix1 = streaksOut->data[3]; // y coordinate for point 1
     350
     351            psArray *chipPix2 = streaksOut->data[4]; // Chip for point 2
     352            psArray *cellPix2 = streaksOut->data[5]; // Cell for point 2
     353            psVector *xPix2 = streaksOut->data[6]; // x coordinate for point 2
     354            psVector *yPix2 = streaksOut->data[7]; // y coordinate for point 2
    326355
    327356            psVector *ra1 = streaks->data[0]; // RA coordinate for point 1
    … …  
    330359            psVector *dec2 = streaks->data[3]; // Dec coordinate for point 2
    331360
    332             psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
    333             psPlane *fp = psPlaneAlloc();    // Focal plane coordinates
    334             psPlane *tp = psPlaneAlloc();    // Tangent plane coordinates
    335             psSphere *sky = psSphereAlloc(); // Sky coordinates
     361            int numCols = psMetadataLookupS32(NULL, hdu->header, "IMNAXIS1"); // Number of columns
     362            int numRows = psMetadataLookupS32(NULL, hdu->header, "IMNAXIS2"); // Number of rows
     363            if (numCols <= 0 || numRows <= 0) {
     364                psError(psErrorCodeLast(), false, "Unable to read size of chip.");
     365                return false;
     366            }
     367
    336368            for (long i = 0; i < xPix1->n; i++) {
    337                 sky->r = ra1->data.F64[i];
    338                 sky->d = dec1->data.F64[i];
    339                 psProject(tp, sky, astromFile->fpa->toSky);
    340                 psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp);
    341                 psPlaneTransformApply(pix, chip->fromFPA, fp);
    342                 xPix1->data.F32[i] = pix->x;
    343                 yPix1->data.F32[i] = pix->y;
    344 
    345                 sky->r = ra2->data.F64[i];
    346                 sky->d = dec2->data.F64[i];
    347                 psProject(tp, sky, astromFile->fpa->toSky);
    348                 psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp);
    349                 psPlaneTransformApply(pix, chip->fromFPA, fp);
    350                 xPix2->data.F32[i] = pix->x;
    351                 yPix2->data.F32[i] = pix->y;
    352             }
    353         }
     369                float x1, y1;   // Coordinates of point 1
     370                coordSky2Chip(&x1, &y1, ra1->data.F64[i], dec1->data.F64[i],
     371                              data->radians, astromFile->fpa, chip);
     372                if ((x1 >= 0 && x1 < numCols && y1 >= 0 && y1 < numRows)) {
     373                    if (rawChip) {
     374                        psString cellName = NULL; // Name of cell
     375                        coordChip2Cell(&cellName, &x1, &y1, x1, y1, cellNames, cellBounds,
     376                                       cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY);
     377                        cellPix1->data[i] = cellName;
     378                    }
     379                    chipPix1->data[i] = psStringCopy(chipName);
     380                    xPix1->data.F32[i] = x1;
     381                    yPix1->data.F32[i] = y1;
     382                }
     383
     384                float x2, y2;   // Coordinates of point 2
     385                coordSky2Chip(&x2, &y2, ra2->data.F64[i], dec2->data.F64[i],
     386                              data->radians, astromFile->fpa, chip);
     387                if ((x2 >= 0 && x2 < numCols && y2 >= 0 && y2 < numRows)) {
     388                    if (rawChip) {
     389                        psString cellName = NULL; // Name of cell
     390                        coordChip2Cell(&cellName, &x2, &y2, x2, y2, cellNames, cellBounds,
     391                                       cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY);
     392                        cellPix2->data[i] = cellName;
     393                    }
     394                    chipPix2->data[i] = psStringCopy(chipName);
     395                    xPix2->data.F32[i] = x2;
     396                    yPix2->data.F32[i] = y2;
     397                }
     398            }
     399        }
     400
     401        psFree(cellNames);
     402        psFree(cellBounds);
     403        psFree(cellParityX);
     404        psFree(cellParityY);
     405        psFree(cellBinX);
     406        psFree(cellBinY);
     407        psFree(cellX0);
     408        psFree(cellY0);
    354409
    355410        // Chip
    … …  
    395450
    396451    if (streaksOut) {
    397         psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1
    398         psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1
    399         psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2
    400         psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2
     452        psArray *chipPix1 = streaksOut->data[0]; // Chip for point 1
     453        psArray *cellPix1 = streaksOut->data[1]; // Cell for point 1
     454        psVector *xPix1 = streaksOut->data[2]; // x coordinate for point 1
     455        psVector *yPix1 = streaksOut->data[3]; // y coordinate for point 1
     456        psArray *chipPix2 = streaksOut->data[4]; // Chip for point 2
     457        psArray *cellPix2 = streaksOut->data[5]; // Cell for point 2
     458        psVector *xPix2 = streaksOut->data[6]; // x coordinate for point 2
     459        psVector *yPix2 = streaksOut->data[7]; // y coordinate for point 2
    401460
    402461        psVector *ra1 = streaks->data[0]; // RA coordinate for point 1
    … …  
    406465
    407466        for (long i = 0; i < xPix1->n; i++) {
    408             if (data->ds9) {
     467            const char *chipName1 = chipPix1->data[i]; // Name of chip for point 1
     468            const char *cellName1 = cellPix1->data[i]; // Name of cell for point 1
     469            const char *chipName2 = chipPix2->data[i]; // Name of chip for point 2
     470            const char *cellName2 = cellPix2->data[i]; // Name of cell for point 2
     471            if (!data->all &&
     472                (!isfinite(xPix1->data.F32[i]) || !isfinite(yPix1->data.F32[i]) ||
     473                 !chipName1 || (rawFile && !cellName1) ||
     474                 !isfinite(xPix2->data.F32[i]) || !isfinite(yPix2->data.F32[i]) ||
     475                 !chipName2 || (rawFile && !cellName2))) {
     476                continue;
     477            }
     478            if (!rawFile && data->ds9) {
     479                // Region file is only appropriate if we're not mapping all the way back to cell coordinates
    409480                fprintf(data->ds9, "image;line(%f,%f,%f,%f) # line= 0 0 color=%s\n",
    410481                        xPix1->data.F32[i], yPix1->data.F32[i], xPix2->data.F32[i], yPix2->data.F32[i],
    411482                        data->ds9color);
    412483            } else {
    413                 fprintf(stdout, "%.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %.3f %.3f\n",
     484                fprintf(stdout, "%.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %s%s%s %.3f %.3f %s%s%s\n",
    414485                        ra1->data.F64[i], dec1->data.F64[i],
    415486                        ra2->data.F64[i], dec2->data.F64[i],
    416487                        xPix1->data.F32[i], yPix1->data.F32[i],
    417                         xPix2->data.F32[i], yPix2->data.F32[i]);
     488                        chipName1 ? chipName1 : "UNKNOWN",
     489                        rawFile ? " " : "",
     490                        rawFile && cellName1 ? cellName1 : (rawFile ? "UNKNOWN" : ""),
     491                        xPix2->data.F32[i], yPix2->data.F32[i],
     492                        chipName2 ? chipName2 : "UNKNOWN",
     493                        rawFile ? " " : "",
     494                        rawFile && cellName2 ? cellName2 : (rawFile ? "UNKNOWN" : "")
     495                    );
    418496            }
    419497        }
Note: See TracChangeset for help on using the changeset viewer.