Index: trunk/ppSkycell/src/ppSkycellCamera.c =================================================================== --- trunk/ppSkycell/src/ppSkycellCamera.c (revision 34089) +++ trunk/ppSkycell/src/ppSkycellCamera.c (revision 34800) @@ -103,11 +103,27 @@ return false; } - if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1")) { + pmFPAfile *bin1 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1"); + if (!bin1) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; } - if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2")) { + pmFPAfile *bin2 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2"); + if (!bin2) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; + } + if (data->masksName) { + pmFPAfile *mask1 = pmFPAfileDefineOutput(data->config, bin1->fpa , "PPSKYCELL.BIN1.MASK"); + if (!mask1) { + psError(psErrorCodeLast(), false, "Unable to define output."); + return false; + } + mask1->save = true; + pmFPAfile *mask2 = pmFPAfileDefineOutput(data->config, bin2->fpa , "PPSKYCELL.BIN2.MASK"); + if (!mask2) { + psError(psErrorCodeLast(), false, "Unable to define output."); + return false; + } + mask2->save = true; } Index: trunk/ppSkycell/src/ppSkycellLoop.c =================================================================== --- trunk/ppSkycell/src/ppSkycellLoop.c (revision 34089) +++ trunk/ppSkycell/src/ppSkycellLoop.c (revision 34800) @@ -15,6 +15,6 @@ ) { - base->x0 = PS_MIN(base->x0, compare->x0); - base->x1 = PS_MAX(base->x1, compare->x1); + base->x0 = PS_MAX(base->x0, compare->x0); + base->x1 = PS_MIN(base->x1, compare->x1); base->y0 = PS_MIN(base->y0, compare->y0); base->y1 = PS_MAX(base->y1, compare->y1); @@ -27,5 +27,5 @@ { psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms - psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell + psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell // Limit the region using the nominated point (X,Y) @@ -38,6 +38,6 @@ toCoords->x += wcs->crpix1; \ toCoords->y += wcs->crpix2; \ - region->x0 = PS_MIN(region->x0, toCoords->x); \ - region->x1 = PS_MAX(region->x1, toCoords->x); \ + region->x0 = PS_MAX(region->x0, toCoords->x); \ + region->x1 = PS_MIN(region->x1, toCoords->x); \ region->y0 = PS_MIN(region->y0, toCoords->y); \ region->y1 = PS_MAX(region->y1, toCoords->y); @@ -132,5 +132,6 @@ psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image - + psArray *regionHDUs = psArrayAlloc(data->numInputs); + // Determine which projection cells we have to deal with. for (int i = 0; i < data->numInputs; i++) { pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine @@ -178,5 +179,5 @@ target->data.S32[i] = numProj; psTrace("ppSkycell", 3, "Image %d uses new projection\n", i); - + psArrayAdd(regionHDUs,1,hdu); numProj++; } @@ -184,10 +185,10 @@ pmFPAfileActivate(data->config->files, false, NULL); - + // Loop over projections for (int i = 0; i < numProj; i++) { psRegion *projRegion = projRegions->data[i]; // Region for skycell projection psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n", i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1); - int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image + int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image // Size of binned image 1 @@ -198,7 +199,7 @@ psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image - psImageInit(image1, 0); - psImageInit(image2, 0); - + psImageInit(image1,NAN); + psImageInit(image2,NAN); + psImage *mask1 = NULL, *mask2 = NULL; // Binned masks if (data->masksName) { @@ -208,10 +209,12 @@ psImageInit(mask2, 0xFF); } - + pmHDU *projhdu = NULL; + int modify_wcs1 = 1; + int modify_wcs2 = 1; + // Loop over inputs to this projection. for (int j = 0; j < data->numInputs; j++) { if (target->data.S32[j] != i) { continue; } - pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j); if (data->masksName) { @@ -227,19 +230,24 @@ pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j); + if (!projhdu) { + projhdu = file->fpa->hdu; + // psMetadataCopy(projhdu->header,file->fpa->hdu->header); + } +#if 0 + else { + if ((psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX1") >= + psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))&& + (psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX2") >= + psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))) { + projhdu = file->fpa->hdu; + psMetadataCopy(projhdu->header,file->fpa->hdu->header); + } + } +#endif + pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input psFree(view); - // Flip images; no idea why this has to be done, but apparently it does - { - psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR); - psFree(inRO->image); - inRO->image = rot; - } - if (inRO->mask) { - psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT); - psFree(inRO->mask); - inRO->mask = rot; - } - + // data->maskVal = 0xffff; pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) { @@ -256,15 +264,14 @@ psRegion *imageRegion = imageRegions->data[j]; // Region for image // Offsets for image on skycell - int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1; - int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1; + int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1; + int yOffset1 = (-imageRegion->y1 + projRegion->y1) / (float)data->bin1; int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; - // XXX Completely neglecting rotations // The skycells are divided up neatly with them all having the same orientation - psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "="); - psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "="); + psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); + psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); if (data->masksName) { - psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "="); - psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "="); + psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M"); + psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M"); } @@ -287,4 +294,8 @@ pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1"); pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2"); + if (data->masksName) { + pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK"); + pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK"); + } } @@ -313,36 +324,161 @@ if (data->doFits) { + + pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0); + // Copy header from projection root hdu + fits1->fpa->hdu = pmHDUAlloc(NULL); + fits1->fpa->hdu->header = psMetadataAlloc(); + psMetadataCopy(fits1->fpa->hdu->header,projhdu->header); + + // Change wcs here +#define WCS_DEBUG 0 + if (modify_wcs1) { + pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); + double cd1f = 1.0 * data->bin1; + double cd2f = 1.0 * data->bin1; +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2 + ); +#endif + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2, + WCS->trans->x->nX, WCS->trans->x->nY, + WCS->trans->y->nX, WCS->trans->y->nY + ); +#endif + for (int q = 0; q <= WCS->trans->x->nX; q++) { + for (int r = 0; r <= WCS->trans->x->nY; r++) { + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); +#if WCS_DEBUG + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->x->coeff[q][r]); +#endif + } + } + for (int q = 0; q <= WCS->trans->y->nX; q++) { + for (int r = 0; r <= WCS->trans->y->nY; r++) { + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); +#if WCS_DEBUG + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->y->coeff[q][r]); +#endif + } + } + pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS); + WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2 + ); +#endif + modify_wcs1 = 0; + } + + + pmChip *Fchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1"); + psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); + psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); + pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1 - pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2 - - // This is a hack to get a functioning header created so the fits images can be written out. +/* // This is a hack to get a functioning header created so the fits images can be written out. */ psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); + + pmReadout *Fro1 = pmReadoutAlloc(Fcell1); + Fro1->image = image1; + Fro1->mask = mask1; + Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; + + fits1->save = true; + fits1->fileIndex = i; + + // Do the mask if we need to + if (data->masksName) { +/* pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */ +/* // Mask1->fpa->hdu = pmHDUAlloc(NULL); */ +/* // Mask1->fpa->hdu->header = psMetadataAlloc(); */ +/* // psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */ +/* pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */ +/* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */ + +/* pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */ +/* /\* // This is a hack to get a functioning header created so the fits images can be written out. *\/ */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */ + +/* pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */ +/* Mro1->image = image1; */ +/* Mro1->mask = mask1; */ +/* Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */ + +/* Mask1->save = true; */ +/* Mask1->fileIndex = i; */ + } + + + + // Repeat with second binned image + pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); + fits2->fpa->hdu = pmHDUAlloc(NULL); + fits2->fpa->hdu->header = psMetadataAlloc(); + psMetadataCopy(fits2->fpa->hdu->header,projhdu->header); + + if (modify_wcs2) { + pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header); + double cd1f = 1.0 * data->bin2 * data->bin1; + double cd2f = 1.0 * data->bin2 * data->bin1; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + for (int q = 0; q < WCS->trans->x->nX; q++) { + for (int r = 0; r < WCS->trans->x->nY; r++) { + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + for (int q = 0; q < WCS->trans->y->nX; q++) { + for (int r = 0; r < WCS->trans->y->nY; r++) { + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + pmAstromWCStoHeader (fits2->fpa->hdu->header,WCS); + modify_wcs2 = 0; + } + + pmChip *Fchip2 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN2"); + psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); + psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); + + + pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); - - psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); - - pmReadout *Fro1 = pmReadoutAlloc(Fcell1), *Fro2 = pmReadoutAlloc(Fcell2); // Binned readouts - - Fro1->image = image1; + + pmReadout *Fro2 = pmReadoutAlloc(Fcell2); Fro2->image = image2; - Fro1->mask = mask1; Fro2->mask = mask2; - - Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true; - - pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0); - fits1->save = true; - fits1->fileIndex = i; - pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); + fits2->save = true; fits2->fileIndex = i; + + }