Index: branches/czw_branch/20120906/ppSkycell/src/ppSkycellLoop.c =================================================================== --- branches/czw_branch/20120906/ppSkycell/src/ppSkycellLoop.c (revision 34410) +++ branches/czw_branch/20120906/ppSkycell/src/ppSkycellLoop.c (revision 34449) @@ -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,5 @@ 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); for (int i = 0; i < data->numInputs; i++) { pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine @@ -178,5 +178,5 @@ target->data.S32[i] = numProj; psTrace("ppSkycell", 3, "Image %d uses new projection\n", i); - + psArrayAdd(regionHDUs,1,hdu); numProj++; } @@ -189,5 +189,5 @@ 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 @@ -208,5 +208,7 @@ psImageInit(mask2, 0xFF); } - + pmHDU *projhdu = NULL; + int modify_wcs1 = 1; + int modify_wcs2 = 1; for (int j = 0; j < data->numInputs; j++) { if (target->data.S32[j] != i) { @@ -227,18 +229,21 @@ pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j); + if (!projhdu) { + projhdu = file->fpa->hdu; + } 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; - } +/* { */ +/* 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; */ +/* } */ pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts @@ -256,5 +261,5 @@ psRegion *imageRegion = imageRegions->data[j]; // Region for image // Offsets for image on skycell - int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1; + int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1; int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1; int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; @@ -313,34 +318,125 @@ 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 + if (modify_wcs1) { + pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); + double cd1f = 1.0 * data->bin1; + double cd2f = 1.0 * data->bin1; + 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 + ); + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + 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 + ); + + 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); + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->x->coeff[q][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); + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->y->coeff[q][r]); + + } + } + pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS); + WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); + 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 + ); + + 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; + + // 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;