Index: trunk/ppViz/src/ppCoord/ppCoordLoop.c =================================================================== --- trunk/ppViz/src/ppCoord/ppCoordLoop.c (revision 27609) +++ trunk/ppViz/src/ppCoord/ppCoordLoop.c (revision 27613) @@ -9,4 +9,61 @@ #include "ppCoord.h" +// Convert sky coordinates to chip coordinates +static void coordSky2Chip(float *x, float *y, // Chip coordinates (output) + double ra, double dec, // Sky coordinates + bool radians, // Coordinates are in radians? + const pmFPA *astromFPA, // Astrometry FPA + const pmChip *astromChip // Astrometry chip + ) +{ + psPlane pix; // Pixel coordinates on chip + psPlane fp; // Focal plane coordinates + psPlane tp; // Tangent plane coordinates + psSphere sky; // Sky coordinates + + sky.r = ra; + sky.d = dec; + + if (!radians) { + sky.r *= M_PI / 180.0; + sky.d *= M_PI / 180.0; + } + + psProject(&tp, &sky, astromFPA->toSky); + psPlaneTransformApply(&fp, astromFPA->fromTPA, &tp); + psPlaneTransformApply(&pix, astromChip->fromFPA, &fp); + + *x = pix.x; + *y = pix.y; + + return; +} + +// Convert chip coordinates to cell coordinates +static void coordChip2Cell(psString *cellName, // Cell name (output) + float *xCell, float *yCell, // Pixel coordinates (output) + float xChip, float yChip, // Sky coordinates + const psArray *cellNames, // Names of cells + const psArray *cellBounds, // Bounds of cells + const psVector *cellX0, const psVector *cellY0, // Cell offsets + const psVector *cellParityX, const psVector *cellParityY, // Cell parities + const psVector *cellBinX, const psVector *cellBinY // Cell binnings + + + ) +{ + int numCells = cellNames->n; // Number of cells + for (int i = 0; i < numCells && !cellName; i++) { + psRegion *region = cellBounds->data[i]; // Bounds of cell + if (xChip >= region->x0 && xChip < region->x1 && yChip >= region->y0 && yChip < region->y1) { + *cellName = psStringCopy(cellNames->data[i]); + // Transform chip coordinates to cell coordinates + *xCell = (xChip - cellX0->data.S32[i]) / (float)(cellParityX->data.S32[i] * cellBinX->data.S32[i]); + *yCell = (yChip - cellY0->data.S32[i]) / (float)(cellParityY->data.S32[i] * cellBinY->data.S32[i]); + } + } + + return; +} bool ppCoordLoop(ppCoordData *data // Run-time data @@ -75,13 +132,17 @@ } } - streaksOut = psArrayAlloc(4); - streaksOut->data[0] = psVectorAlloc(numStreaks, PS_TYPE_F32); - streaksOut->data[1] = psVectorAlloc(numStreaks, PS_TYPE_F32); + streaksOut = psArrayAlloc(8); + streaksOut->data[0] = psArrayAlloc(numStreaks); + streaksOut->data[1] = psArrayAlloc(numStreaks); streaksOut->data[2] = psVectorAlloc(numStreaks, PS_TYPE_F32); streaksOut->data[3] = psVectorAlloc(numStreaks, PS_TYPE_F32); - psVectorInit(streaksOut->data[0], NAN); - psVectorInit(streaksOut->data[1], NAN); - psVectorInit(streaksOut->data[2], NAN); + streaksOut->data[4] = psArrayAlloc(numStreaks); + streaksOut->data[5] = psArrayAlloc(numStreaks); + streaksOut->data[6] = psVectorAlloc(numStreaks, PS_TYPE_F32); + streaksOut->data[7] = psVectorAlloc(numStreaks, PS_TYPE_F32); psVectorInit(streaksOut->data[3], NAN); + psVectorInit(streaksOut->data[4], NAN); + psVectorInit(streaksOut->data[6], NAN); + psVectorInit(streaksOut->data[7], NAN); } @@ -183,66 +244,66 @@ } + + pmChip *rawChip = rawFile ? pmFPAviewThisChip(view, rawFile->fpa) : NULL; // Chip with raw + psArray *cellBounds = NULL; // Bounds of each cell + psArray *cellNames = NULL; // Names of each cell + psVector *cellParityX = NULL, *cellParityY = NULL; // Parity for each cell + psVector *cellX0 = NULL, *cellY0 = NULL; // Offset for each cell + psVector *cellBinX = NULL, *cellBinY = NULL; // Binning for each cell + if (rawChip) { + if (!rawChip->data_exists) { + // Not interested in this chip + continue; + } + int numCells = rawChip->cells->n; // Number of cells + + cellBounds = psArrayAlloc(numCells); + cellNames = psArrayAlloc(numCells); + cellParityX = psVectorAlloc(numCells, PS_TYPE_S32); + cellParityY = psVectorAlloc(numCells, PS_TYPE_S32); + cellX0 = psVectorAlloc(numCells, PS_TYPE_S32); + cellY0 = psVectorAlloc(numCells, PS_TYPE_S32); + cellBinX = psVectorAlloc(numCells, PS_TYPE_S32); + cellBinY = psVectorAlloc(numCells, PS_TYPE_S32); + + pmCell *rawCell; // Cell with raw data + while ((rawCell = pmFPAviewNextCell(view, rawFile->fpa, 1))) { + if (!rawCell->data_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Error loading data from files."); + return false; + } + psRegion *region = pmCellExtent(rawCell); // Bounds of cell + if (!region) { + psError(psErrorCodeLast(), false, "Unable to determine extent of cell."); + return false; + } + cellBounds->data[view->cell] = region; + cellNames->data[view->cell] = psMemIncrRefCounter(psMetadataLookupStr(NULL, rawCell->concepts, "CELL.NAME")); + cellParityX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XPARITY"); + cellParityY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YPARITY"); + cellX0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.X0"); + cellY0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.Y0"); + cellBinX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XBIN"); + cellBinY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YBIN"); + if (cellParityX->data.S32[view->cell] == 0 || cellParityY->data.S32[view->cell] == 0 || + cellBinX->data.S32[view->cell] == 0 || cellBinY->data.S32[view->cell] == 0) { + psError(PM_ERR_CONCEPTS, true, "Concepts aren't set: %d %d %d %d %d %d\n", + cellParityX->data.S32[view->cell], cellParityY->data.S32[view->cell], + cellX0->data.S32[view->cell], cellY0->data.S32[view->cell], + cellBinX->data.S32[view->cell], cellBinY->data.S32[view->cell]); + return false; + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Error freeing data from files."); + return false; + } + } + } + if (radec) { - pmChip *rawChip = rawFile ? pmFPAviewThisChip(view, rawFile->fpa) : NULL; // Chip with raw - psArray *cellBounds = NULL; // Bounds of each cell - psArray *cellNames = NULL; // Names of each cell - psVector *cellParityX = NULL, *cellParityY = NULL; // Parity for each cell - psVector *cellX0 = NULL, *cellY0 = NULL; // Offset for each cell - psVector *cellBinX = NULL, *cellBinY = NULL; // Binning for each cell - if (rawChip) { - if (!rawChip->data_exists) { - // Not interested in this chip - continue; - } - int numCells = rawChip->cells->n; // Number of cells - - cellBounds = psArrayAlloc(numCells); - cellNames = psArrayAlloc(numCells); - cellParityX = psVectorAlloc(numCells, PS_TYPE_S32); - cellParityY = psVectorAlloc(numCells, PS_TYPE_S32); - cellX0 = psVectorAlloc(numCells, PS_TYPE_S32); - cellY0 = psVectorAlloc(numCells, PS_TYPE_S32); - cellBinX = psVectorAlloc(numCells, PS_TYPE_S32); - cellBinY = psVectorAlloc(numCells, PS_TYPE_S32); - - pmCell *rawCell; // Cell with raw data - while ((rawCell = pmFPAviewNextCell(view, rawFile->fpa, 1))) { - if (!rawCell->data_exists) { - continue; - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(psErrorCodeLast(), false, "Error loading data from files."); - return false; - } - psRegion *region = pmCellExtent(rawCell); // Bounds of cell - if (!region) { - psError(psErrorCodeLast(), false, "Unable to determine extent of cell."); - return false; - } - cellBounds->data[view->cell] = region; - cellNames->data[view->cell] = psMemIncrRefCounter(psMetadataLookupStr(NULL, - rawCell->concepts, "CELL.NAME")); - cellParityX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XPARITY"); - cellParityY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YPARITY"); - cellX0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.X0"); - cellY0->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.Y0"); - cellBinX->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.XBIN"); - cellBinY->data.S32[view->cell] = psMetadataLookupS32(NULL, rawCell->concepts, "CELL.YBIN"); - if (cellParityX->data.S32[view->cell] == 0 || cellParityY->data.S32[view->cell] == 0 || - cellBinX->data.S32[view->cell] == 0 || cellBinY->data.S32[view->cell] == 0) { - psError(PM_ERR_CONCEPTS, true, "Concepts aren't set: %d %d %d %d %d %d\n", - cellParityX->data.S32[view->cell], cellParityY->data.S32[view->cell], - cellX0->data.S32[view->cell], cellY0->data.S32[view->cell], - cellBinX->data.S32[view->cell], cellBinY->data.S32[view->cell]); - return false; - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "Error freeing data from files."); - return false; - } - } - } - psVector *ra = radec->data[0], *dec = radec->data[1]; // Pixel coordinates long num = ra->n; // Number of coordinates @@ -255,9 +316,4 @@ } - psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip - psPlane *fp = psPlaneAlloc(); // Focal plane coordinates - psPlane *tp = psPlaneAlloc(); // Tangent plane coordinates - psSphere *sky = psSphereAlloc(); // Sky coordinates - psArray *chipPix = radecOut->data[0]; // Chip for pixels psVector *xPix = radecOut->data[1]; // x coordinate for pixels @@ -266,17 +322,7 @@ for (long i = 0; i < num; i++) { - sky->r = ra->data.F64[i]; - sky->d = dec->data.F64[i]; - - if (!data->radians) { - sky->r *= M_PI / 180.0; - sky->d *= M_PI / 180.0; - } - - psProject(tp, sky, astromFile->fpa->toSky); - psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp); - psPlaneTransformApply(pix, chip->fromFPA, fp); - - float x = pix->x, y = pix->y; // Pixel coordinates + float x, y; // Pixel coordinates + coordSky2Chip(&x, &y, ra->data.F64[i], dec->data.F64[i], + data->radians, astromFile->fpa, chip); if ((x < 0 || x > numCols || y < 0 || y > numRows)) { // Not on this chip @@ -285,17 +331,7 @@ if (rawChip) { - psString cellName = NULL; // Name of cell - int numCells = rawChip->cells->n; // Number of cells - for (int i = 0; i < numCells && !cellName; i++) { - psRegion *region = cellBounds->data[i]; // Bounds of cell - if (x >= region->x0 && x < region->x1 && y >= region->y0 && y < region->y1) { - cellName = psStringCopy(cellNames->data[i]); - // Transform chip coordinates to cell coordinates - x = (x - cellX0->data.S32[i]) / - (float)(cellParityX->data.S32[i] * cellBinX->data.S32[i]); - y = (y - cellY0->data.S32[i]) / - (float)(cellParityY->data.S32[i] * cellBinY->data.S32[i]); - } - } + psString cellName = NULL; // Name of cell + coordChip2Cell(&cellName, &x, &y, x, y, cellNames, cellBounds, + cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY); cellPix->data[i] = cellName; } @@ -305,23 +341,16 @@ yPix->data.F32[i] = y; } - psFree(pix); - psFree(fp); - psFree(tp); - psFree(sky); - psFree(cellNames); - psFree(cellBounds); - psFree(cellParityX); - psFree(cellParityY); - psFree(cellBinX); - psFree(cellBinY); - psFree(cellX0); - psFree(cellY0); - } - - if (streaks && data->chipName && !rawFile) { - psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1 - psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1 - psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2 - psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2 + } + + if (streaks) { + psArray *chipPix1 = streaksOut->data[0]; // Chip for point 1 + psArray *cellPix1 = streaksOut->data[1]; // Cell for point 1 + psVector *xPix1 = streaksOut->data[2]; // x coordinate for point 1 + psVector *yPix1 = streaksOut->data[3]; // y coordinate for point 1 + + psArray *chipPix2 = streaksOut->data[4]; // Chip for point 2 + psArray *cellPix2 = streaksOut->data[5]; // Cell for point 2 + psVector *xPix2 = streaksOut->data[6]; // x coordinate for point 2 + psVector *yPix2 = streaksOut->data[7]; // y coordinate for point 2 psVector *ra1 = streaks->data[0]; // RA coordinate for point 1 @@ -330,26 +359,52 @@ psVector *dec2 = streaks->data[3]; // Dec coordinate for point 2 - psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip - psPlane *fp = psPlaneAlloc(); // Focal plane coordinates - psPlane *tp = psPlaneAlloc(); // Tangent plane coordinates - psSphere *sky = psSphereAlloc(); // Sky 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; + } + for (long i = 0; i < xPix1->n; i++) { - sky->r = ra1->data.F64[i]; - sky->d = dec1->data.F64[i]; - psProject(tp, sky, astromFile->fpa->toSky); - psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp); - psPlaneTransformApply(pix, chip->fromFPA, fp); - xPix1->data.F32[i] = pix->x; - yPix1->data.F32[i] = pix->y; - - sky->r = ra2->data.F64[i]; - sky->d = dec2->data.F64[i]; - psProject(tp, sky, astromFile->fpa->toSky); - psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp); - psPlaneTransformApply(pix, chip->fromFPA, fp); - xPix2->data.F32[i] = pix->x; - yPix2->data.F32[i] = pix->y; - } - } + float x1, y1; // Coordinates of point 1 + coordSky2Chip(&x1, &y1, ra1->data.F64[i], dec1->data.F64[i], + data->radians, astromFile->fpa, chip); + if ((x1 >= 0 && x1 < numCols && y1 >= 0 && y1 < numRows)) { + if (rawChip) { + psString cellName = NULL; // Name of cell + coordChip2Cell(&cellName, &x1, &y1, x1, y1, cellNames, cellBounds, + cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY); + cellPix1->data[i] = cellName; + } + chipPix1->data[i] = psStringCopy(chipName); + xPix1->data.F32[i] = x1; + yPix1->data.F32[i] = y1; + } + + float x2, y2; // Coordinates of point 2 + coordSky2Chip(&x2, &y2, ra2->data.F64[i], dec2->data.F64[i], + data->radians, astromFile->fpa, chip); + if ((x2 >= 0 && x2 < numCols && y2 >= 0 && y2 < numRows)) { + if (rawChip) { + psString cellName = NULL; // Name of cell + coordChip2Cell(&cellName, &x2, &y2, x2, y2, cellNames, cellBounds, + cellX0, cellY0, cellParityX, cellParityY, cellBinX, cellBinY); + cellPix2->data[i] = cellName; + } + chipPix2->data[i] = psStringCopy(chipName); + xPix2->data.F32[i] = x2; + yPix2->data.F32[i] = y2; + } + } + } + + psFree(cellNames); + psFree(cellBounds); + psFree(cellParityX); + psFree(cellParityY); + psFree(cellBinX); + psFree(cellBinY); + psFree(cellX0); + psFree(cellY0); // Chip @@ -395,8 +450,12 @@ if (streaksOut) { - psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1 - psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1 - psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2 - psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2 + psArray *chipPix1 = streaksOut->data[0]; // Chip for point 1 + psArray *cellPix1 = streaksOut->data[1]; // Cell for point 1 + psVector *xPix1 = streaksOut->data[2]; // x coordinate for point 1 + psVector *yPix1 = streaksOut->data[3]; // y coordinate for point 1 + psArray *chipPix2 = streaksOut->data[4]; // Chip for point 2 + psArray *cellPix2 = streaksOut->data[5]; // Cell for point 2 + psVector *xPix2 = streaksOut->data[6]; // x coordinate for point 2 + psVector *yPix2 = streaksOut->data[7]; // y coordinate for point 2 psVector *ra1 = streaks->data[0]; // RA coordinate for point 1 @@ -406,14 +465,33 @@ for (long i = 0; i < xPix1->n; i++) { - if (data->ds9) { + const char *chipName1 = chipPix1->data[i]; // Name of chip for point 1 + const char *cellName1 = cellPix1->data[i]; // Name of cell for point 1 + const char *chipName2 = chipPix2->data[i]; // Name of chip for point 2 + const char *cellName2 = cellPix2->data[i]; // Name of cell for point 2 + if (!data->all && + (!isfinite(xPix1->data.F32[i]) || !isfinite(yPix1->data.F32[i]) || + !chipName1 || (rawFile && !cellName1) || + !isfinite(xPix2->data.F32[i]) || !isfinite(yPix2->data.F32[i]) || + !chipName2 || (rawFile && !cellName2))) { + 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;line(%f,%f,%f,%f) # line= 0 0 color=%s\n", xPix1->data.F32[i], yPix1->data.F32[i], xPix2->data.F32[i], yPix2->data.F32[i], data->ds9color); } else { - fprintf(stdout, "%.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %.3f %.3f\n", + fprintf(stdout, "%.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], xPix1->data.F32[i], yPix1->data.F32[i], - xPix2->data.F32[i], yPix2->data.F32[i]); + chipName1 ? chipName1 : "UNKNOWN", + rawFile ? " " : "", + rawFile && cellName1 ? cellName1 : (rawFile ? "UNKNOWN" : ""), + xPix2->data.F32[i], yPix2->data.F32[i], + chipName2 ? chipName2 : "UNKNOWN", + rawFile ? " " : "", + rawFile && cellName2 ? cellName2 : (rawFile ? "UNKNOWN" : "") + ); } }