IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 4, 2021, 6:14:46 PM (5 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-dev-20210817 (speed-ups for inverse transform, option for extra orders)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r32600 r41899  
    99#include "ppCoord.h"
    1010
     11typedef struct {
     12  psPlane pix;
     13  psPlane fp;
     14  psPlane tp;
     15  psSphere sky;
     16} ppCoordSet;
     17
    1118// Convert sky coordinates to chip coordinates
    12 static void coordSky2Chip(float *x, float *y, // Chip coordinates (output)
    13                           double ra, double dec, // Sky coordinates
     19static void coordSky2Chip(ppCoordSet *src, // Sky coordinates
    1420                          bool radians,          // Coordinates are in radians?
    1521                          const pmFPA *astromFPA,  // Astrometry FPA
    … …  
    1723    )
    1824{
    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 
     25    // convert from src->sky to src->pix
    2726    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     psTrace("ppCoord",2,"sky2pix tr1: %f %f\n",tp.x,tp.y);
    34     psPlaneTransformApply(&fp, astromFPA->fromTPA, &tp);
    35     psTrace("ppCoord",2,"sky2pix tr2: %f %f\n",fp.x,fp.y);
    36     psPlaneTransformApply(&pix, astromChip->fromFPA, &fp);
    37     psTrace("ppCoord",2,"sky2pix tr3: %f %f\n",pix.x,pix.y);
    38 
    39     *x = pix.x;
    40     *y = pix.y;
     27        src->sky.r *= M_PI / 180.0;
     28        src->sky.d *= M_PI / 180.0;
     29    }
     30
     31    psProject(&src->tp, &src->sky, astromFPA->toSky);
     32    psTrace("ppCoord",2,"sky2pix tr1: %f %f\n", src->tp.x, src->tp.y);
     33    psPlaneTransformApply(&src->fp, astromFPA->fromTPA, &src->tp);
     34    psTrace("ppCoord",2,"sky2pix tr2: %f %f\n", src->fp.x, src->fp.y);
     35    psPlaneTransformApply(&src->pix, astromChip->fromFPA, &src->fp);
     36    psTrace("ppCoord",2,"sky2pix tr3: %f %f\n", src->pix.x, src->pix.y);
    4137
    4238    return;
    … …  
    10096        psVector *ra = radec->data[0];  // RA coordinates
    10197        long num = ra->n;               // Number of coordinates
    102         radecOut = psArrayAlloc(4);
     98        radecOut = psArrayAlloc(8);
    10399        radecOut->data[0] = psArrayAlloc(num);
    104100        radecOut->data[1] = psVectorAlloc(num, PS_TYPE_F32);
    105101        radecOut->data[2] = psVectorAlloc(num, PS_TYPE_F32);
    106102        radecOut->data[3] = psArrayAlloc(num);
     103        radecOut->data[4] = psVectorAlloc(num, PS_TYPE_F32); // FP x
     104        radecOut->data[5] = psVectorAlloc(num, PS_TYPE_F32); // FP y
     105        radecOut->data[6] = psVectorAlloc(num, PS_TYPE_F32); // TP x
     106        radecOut->data[7] = psVectorAlloc(num, PS_TYPE_F32); // TP y
    107107        psVectorInit(radecOut->data[1], NAN);
    108108        psVectorInit(radecOut->data[2], NAN);
    … …  
    229229        }
    230230
     231        // fprintf (stderr, "chip: %s\n", chipName);
     232
    231233        // read WCS data from the corresponding header
    232234        pmHDU *hdu = pmFPAviewThisHDU (view, astromFile->fpa);
    … …  
    274276                }
    275277
    276                 fprintf(stdout, "%s %.3f %.3f --> %.10lf %.10lf\n", chipName, pix->x, pix->y, sky->r, sky->d);
     278                fprintf(stdout, "%s %10.3f %10.3f --> %14.10lf %14.10lf : %11.4f %11.4f : %11.4f %11.4f\n", chipName, pix->x, pix->y, sky->r, sky->d, tp->x, tp->y, fp->x, fp->y);
    277279            }
    278280            psFree(pix);
    … …  
    358360            psVector *yPix = radecOut->data[2];   // y coordinate for pixels
    359361            psArray *cellPix = radecOut->data[3]; // Cell for pixels
     362            psVector *xFP = radecOut->data[4];   // x coordinate for FP
     363            psVector *yFP = radecOut->data[5];   // y coordinate for FP
     364            psVector *xTP = radecOut->data[6];   // x coordinate for TP
     365            psVector *yTP = radecOut->data[7];   // y coordinate for TP
    360366
    361367            for (long i = 0; i < num; i++) {
    362                 float x, y;             // Pixel coordinates
    363                 coordSky2Chip(&x, &y, ra->data.F64[i], dec->data.F64[i],
    364                                 data->radians, astromFile->fpa, chip);
    365                 if ((x < 0 || x > numCols || y < 0 || y > numRows)) {
     368                ppCoordSet src;
     369                src.sky.r = ra->data.F64[i];
     370                src.sky.d = dec->data.F64[i];
     371                coordSky2Chip(&src, data->radians, astromFile->fpa, chip);
     372                // fprintf (stderr, "x,y: %.2f %.2f\n", src.pix.x, src.pix.y);
     373                if ((src.pix.x < 0 || src.pix.x > numCols || src.pix.y < 0 || src.pix.y > numRows)) {
    366374                    // Not on this chip
    367375                    continue;
    … …  
    369377
    370378                if (rawChip) {
     379                    float x, y;
    371380                    psString cellName = NULL; // Name of cell
    372                     coordChip2Cell(&cellName, &x, &y, x, y, cellNames, cellBounds,
     381                    coordChip2Cell(&cellName, &x, &y, src.pix.x, src.pix.y, cellNames, cellBounds,
    373382                                   cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY);
    374383                    cellPix->data[i] = cellName;
    375                 }
    376 
     384                    xPix->data.F32[i] = x;
     385                    yPix->data.F32[i] = y;
     386                } else {
     387                    xPix->data.F32[i] = src.pix.x;
     388                    yPix->data.F32[i] = src.pix.y;
     389                }
    377390                chipPix->data[i] = psStringCopy(chipName);
    378                 xPix->data.F32[i] = x;
    379                 yPix->data.F32[i] = y;
     391                xFP->data.F32[i] = src.fp.x;
     392                yFP->data.F32[i] = src.fp.y;
     393                xTP->data.F32[i] = src.tp.x;
     394                yTP->data.F32[i] = src.tp.y;
    380395            }
    381396        }
    … …  
    406421            for (long i = 0; i < xPix1->n; i++) {
    407422                float x1, y1;   // Coordinates of point 1
    408                 coordSky2Chip(&x1, &y1, ra1->data.F64[i], dec1->data.F64[i],
    409                               true, astromFile->fpa, chip);
     423                ppCoordSet src;
     424                src.sky.r = ra1->data.F64[i];
     425                src.sky.d = dec1->data.F64[i];
     426                coordSky2Chip(&src, true, astromFile->fpa, chip);
     427                x1 = src.pix.x;
     428                y1 = src.pix.y;
    410429                if ((x1 >= 0 && x1 < numCols && y1 >= 0 && y1 < numRows)) {
    411430                    if (rawChip) {
    … …  
    421440
    422441                float x2, y2;   // Coordinates of point 2
    423                 coordSky2Chip(&x2, &y2, ra2->data.F64[i], dec2->data.F64[i],
    424                               true, astromFile->fpa, chip);
     442                src.sky.r = ra2->data.F64[i];
     443                src.sky.d = dec2->data.F64[i];
     444                coordSky2Chip(&src, true, astromFile->fpa, chip);
     445                x2 = src.pix.x;
     446                y2 = src.pix.y;
    425447                if ((x2 >= 0 && x2 < numCols && y2 >= 0 && y2 < numRows)) {
    426448                    if (rawChip) {
    … …  
    455477            for (long i = 0; i < num; i++) {
    456478                float x, y;             // Pixel coordinates
    457                 coordSky2Chip(&x, &y, ra->data.F64[i], dec->data.F64[i],
    458                                 false, astromFile->fpa, chip);
     479                ppCoordSet src;
     480                src.sky.r = ra->data.F64[i];
     481                src.sky.d = dec->data.F64[i];
     482                coordSky2Chip(&src, false, astromFile->fpa, chip);
     483                x = src.pix.x;
     484                y = src.pix.y;
    459485                if ((x < 0 || x > numCols || y < 0 || y > numRows)) {
    460486                    // Not on this chip
    … …  
    503529        psVector *ra = radec->data[0];        // RA coordinate
    504530        psVector *dec = radec->data[1];       // Dec coordinate
     531        psVector *xFP = radecOut->data[4];   // x coordinate for FP
     532        psVector *yFP = radecOut->data[5];   // y coordinate for FP
     533        psVector *xTP = radecOut->data[6];   // x coordinate for TP
     534        psVector *yTP = radecOut->data[7];   // y coordinate for TP
    505535
    506536        for (long i = 0; i < chipPix->n; i++) {
    … …  
    516546                        xPix->data.F32[i], yPix->data.F32[i], data->ds9radius, data->ds9color);
    517547            } else {
    518                 fprintf(stdout, "%.10lf %.10lf --> %.3f %.3f %s%s%s\n",
     548                fprintf(stdout, "%14.10lf %14.10lf --> %10.3f %10.3f %s%s%s : %11.4f %11.4f : %11.4f %11.4f\n",
    519549                        ra->data.F64[i], dec->data.F64[i], xPix->data.F32[i], yPix->data.F32[i],
    520550                        chipName ? chipName : "UNKNOWN",
    521551                        rawFile ? " " : "",
    522                         rawFile && cellName ? cellName : (rawFile ? "UNKNOWN" : ""));
     552                        rawFile && cellName ? cellName : (rawFile ? "UNKNOWN" : ""),
     553                        xFP->data.F32[i], yFP->data.F32[i], xTP->data.F32[i], yTP->data.F32[i] );
    523554            }
    524555        }
Note: See TracChangeset for help on using the changeset viewer.