Index: /trunk/ppSkycell/src/ppSkycell.h =================================================================== --- /trunk/ppSkycell/src/ppSkycell.h (revision 36356) +++ /trunk/ppSkycell/src/ppSkycell.h (revision 36357) @@ -10,5 +10,4 @@ typedef struct { psString imagesName; // Filename with images - psString masksName; // Filename with masks psString wcsrefName; // Filename with WCS references psString outRoot; // Output root name Index: /trunk/ppSkycell/src/ppSkycellArguments.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellArguments.c (revision 36356) +++ /trunk/ppSkycell/src/ppSkycellArguments.c (revision 36357) @@ -29,5 +29,5 @@ { fprintf(stderr, "\nPan-STARRS skycell JPEGifier\n\n"); - fprintf(stderr, "Usage: %s -images INPUT.list [-masks MASK.list] OUTPUT_ROOT\n\n", + fprintf(stderr, "Usage: %s -images INPUT.list OUTPUT_ROOT\n\n", program); fprintf(stderr, "\n"); @@ -50,5 +50,4 @@ psMetadata *arguments = psMetadataAlloc(); // Command-line arguments psMetadataAddStr(arguments, PS_LIST_TAIL, "-images", 0, "Filename with input images", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-masks", 0, "Filename with input masks", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-wcsref", 0, "Filename with WCS references", NULL); psMetadataAddS16(arguments, PS_LIST_TAIL, "-exptimeOrder", 0, "Order of exptime scaling", 0); @@ -58,5 +57,4 @@ data->imagesName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-images")); - data->masksName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-masks")); data->wcsrefName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-wcsref")); data->exptimeOrder = psMetadataLookupS16(NULL,arguments,"-exptimeOrder"); Index: /trunk/ppSkycell/src/ppSkycellCamera.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellCamera.c (revision 36356) +++ /trunk/ppSkycell/src/ppSkycellCamera.c (revision 36357) @@ -56,21 +56,4 @@ data->numInputs = images->n; - psArray *masks = NULL; // Mask names - if (data->masksName) { - masks = fileList(data->masksName); - if (!masks) { - psError(psErrorCodeLast(), false, "No masks provided."); - psFree(images); - return false; - } - if (masks->n != data->numInputs) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of images (%ld) and masks (%ld) do not match", - images->n, masks->n); - psFree(images); - psFree(masks); - return false; - } - } - psArray *wcsref = NULL; // Names of WCS reference images if (data->wcsrefName) { @@ -79,5 +62,4 @@ psError(psErrorCodeLast(), false, "No WCSrefs provided."); psFree(images); - psFree(masks); return(false); } @@ -86,5 +68,4 @@ images->n, wcsref->n); psFree(images); - psFree(masks); psFree(wcsref); return(false); @@ -104,13 +85,4 @@ } - if (data->masksName) { - fileArguments("MASK", masks->data[i], "Name of the mask", data->config); - if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.MASK", "MASK") || !status) { - psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.MASK"); - // XXX Cleanup - return false; - } - } - if (data->wcsrefName) { fileArguments("WCSREF", wcsref->data[i], "Name of the WCS reference", data->config); @@ -120,10 +92,5 @@ return(false); } -/* if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.WCSREF", "WCSREF") || !status) { */ -/* psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.WCSREF"); */ -/* return false; */ -/* } */ } - } @@ -145,18 +112,4 @@ 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/ppSkycellData.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellData.c (revision 36356) +++ /trunk/ppSkycell/src/ppSkycellData.c (revision 36357) @@ -14,5 +14,4 @@ { psFree(data->imagesName); - psFree(data->masksName); psFree(data->wcsrefName); psFree(data->outRoot); @@ -28,5 +27,4 @@ data->imagesName = NULL; - data->masksName = NULL; data->wcsrefName = NULL; data->outRoot = NULL; Index: /trunk/ppSkycell/src/ppSkycellLoop.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellLoop.c (revision 36356) +++ /trunk/ppSkycell/src/ppSkycellLoop.c (revision 36357) @@ -11,51 +11,4 @@ #define BUFFER 16 // Size of buffer for projections -#if (0) -// This is all junk that is basically useless, because it calculates things we already know. - -static void regionMinMax(psRegion *base,// Base region; modified - const psRegion *compare // Comparison region - ) -{ - 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); - return; -} - - -static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System - int numCols, int numRows // Size of image - ) -{ - psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms - psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell - -// Limit the region using the nominated point (X,Y) -#define REGION_RANGE(X,Y) \ - fromCoords->x = (X); \ - fromCoords->y = (Y); \ - psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \ - printf("%f %f %f %f\n",toCoords->x,toCoords->y,fromCoords->x,fromCoords->y); \ - toCoords->x /= wcs->cdelt1; \ - toCoords->y /= wcs->cdelt2; \ - toCoords->x -= wcs->crpix1; \ - toCoords->y += wcs->crpix2; \ - 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); - // - REGION_RANGE(0, 0); // Lower left - psTrace("ppSkycell", 9, "Start: %.0f %.0f", toCoords->x, toCoords->y); - REGION_RANGE(0, numRows); // Upper left - REGION_RANGE(numCols, 0); // Lower right - REGION_RANGE(numCols, numRows); // Upper right - psTrace("ppSkycell", 9, "Stop: %.0f %.0f\n", toCoords->x, toCoords->y); - psTrace("ppSkycell", 7, "%g %g %g %g\n",wcs->cdelt1,wcs->cdelt2,wcs->crpix1,wcs->crpix2); - return region; -} -#endif // Activate/deactivate a single element for a list @@ -140,5 +93,5 @@ psArray *regionHDUs = psArrayAlloc(data->numInputs); - psVector *exptimes = psVectorAllocEmpty(data->numInputs, PS_TYPE_F64); + psVector *exptimes = psVectorAlloc(data->numInputs, PS_TYPE_F64); // Determine which projection cells we have to deal with. @@ -146,10 +99,4 @@ pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine // Header in the FPA should have been read as a part of defining the file... -#if 0 - bool pmFPAReadHeaderSet(pmFPA *fpa, // FPA to read into - psFits *fits, // FITS file from which to read - pmConfig *config // Configuration - ); -#endif pmHDU *hdu = file->fpa->hdu; // Header of interest @@ -218,29 +165,27 @@ for (int i = 0; i < numProj; i++) { psRegion *projRegion = projRegions->data[i]; // Region for skycell projection - // projRegion->x0 += 5760; - //projRegion->x1 += 5760; - // projRegion->x0 += - // projRegion->x1 += 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 ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image - // Size of binned image 1 + + // Size of unbinned image + int xSize = projRegion->x1 - projRegion->x0 + 1; + int ySize = projRegion->y1 - projRegion->y0 + 1; + + // Size of binned images int numCols1 = xSize / (float)data->bin1 + 1.5, numRows1 = ySize / (float)data->bin1 + 1.5; - // Size of binned image 2 int numCols2 = numCols1 / (float)data->bin2 + 1.5, numRows2 = numRows1 / (float)data->bin2 + 1.5; - psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image - psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image + // Binned image containers + psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); + psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); psImageInit(image1,NAN); psImageInit(image2,NAN); + + // Binned image radius values. Used to determine primacy. + psImage *radius1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); + psImage *radius2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); + psImageInit(radius1,pow(numCols1 + numRows1,2)); // These values can be anything, just need to be larger than a binned radius. + psImageInit(radius2,pow(numCols1 + numRows1,2)); - psImage *mask1 = NULL, *mask2 = NULL; // Binned masks - if (data->masksName) { - mask1 = psImageAlloc(numCols1, numRows1, PS_TYPE_IMAGE_MASK); - mask2 = psImageAlloc(numCols2, numRows2, PS_TYPE_IMAGE_MASK); - psImageInit(mask1, 0xFF); - psImageInit(mask2, 0xFF); - } // HDU containing the WCS we plan on using pmHDU *projhdu = NULL; @@ -260,7 +205,4 @@ } pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j); - if (data->masksName) { - pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.MASK", j); - } pmFPAview *view = filesIterateDown(data->config); // View to readout @@ -317,40 +259,73 @@ // Scale by exposure time double exptime_factor = pow(exptimes->data.F64[j]/stats->robustMedian,data->exptimeOrder); - psBinaryOp(bin1RO->image,bin1RO->image,"/", - psScalarAlloc(exptime_factor, PS_TYPE_F64)); - psBinaryOp(bin2RO->image,bin2RO->image,"/", - psScalarAlloc(exptime_factor, PS_TYPE_F32)); + if ((isfinite(exptime_factor))&&(exptime_factor != 0.0)) { + psBinaryOp(bin1RO->image,bin1RO->image,"/", + psScalarAlloc(exptime_factor, PS_TYPE_F32)); + psBinaryOp(bin2RO->image,bin2RO->image,"/", + psScalarAlloc(exptime_factor, PS_TYPE_F32)); + } // End scaling -#if (0) - psRegion *imageRegion = imageRegions->data[j]; // Region for image - fprintf(stderr,"%d IM %f %f %f %f PR %f %f %f %f\n", - j, - imageRegion->x0,imageRegion->x1, - imageRegion->y0,imageRegion->y1, - projRegion->x0,projRegion->x1, - projRegion->y0,projRegion->y1); -#endif + // Offsets for image on this projection cell are just differences in CRPIX positions. int xOffset1 = (-1 * wcs->crpix1 - projRegion->x0) / (float)(data->bin1); int yOffset1 = (-1 * wcs->crpix2 - projRegion->y0) / (float)(data->bin1); - int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; + int xOffset2 = xOffset1 / (float)data->bin2; + int yOffset2 = yOffset1 / (float)data->bin2; psTrace("ppSkycell",5,"Offsets: %d %d : %d %d", xOffset1,yOffset1,xOffset2,yOffset2); - // This is where any primary code would go. I'm unclear on how it would work. - + // Check each pixel for primacy. A pixel is primary + // if it is closer to the central pixel of its skycell + // than any other pixel is to theirs. + + // Let's just do the overlay here, instead of doing math, then handing it + // off to another function. That seems silly. + + int u,v,x,y; + + double u0 = xOffset1 + (bin1RO->image->numCols) / 2.0; + double v0 = yOffset1 + (bin1RO->image->numRows) / 2.0; + for (x = 0; x < bin1RO->image->numCols; x++) { + for (y = 0; y < bin1RO->image->numRows; y++) { + if (!isfinite(bin1RO->image->data.F32[y][x])) { + continue; + } + u = x + xOffset1; + v = y + yOffset1; + double R2 = pow(u - u0,2) + pow(v - v0,2); + + if (R2 < radius1->data.F32[v][u]) { + radius1->data.F32[v][u] = R2; + image1->data.F32[v][u] = bin1RO->image->data.F32[y][x]; + } + } + } + + u0 = xOffset2 + (bin2RO->image->numCols) / 2.0; + v0 = yOffset2 + (bin2RO->image->numRows) / 2.0; + for (x = 0; x < bin2RO->image->numCols; x++) { + for (y = 0; y < bin2RO->image->numRows; y++) { + if (!isfinite(bin2RO->image->data.F32[y][x])) { + continue; + } + u = x + xOffset2; + v = y + yOffset2; + double R2 = pow(u - u0,2) + pow(v - v0,2); + if (R2 < radius2->data.F32[v][u]) { + radius2->data.F32[v][u] = R2; + image2->data.F32[v][u] = bin2RO->image->data.F32[y][x]; + } + } + } - + + // Overlay the data onto the appropriate pixels in the final outputs // XXX Completely neglecting rotations // The skycells are divided up neatly with them all having the same orientation - psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); - psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); - if (data->masksName) { - psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M"); - psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M"); - } + // psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); + // psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); // Cleanup on input loop. @@ -360,9 +335,4 @@ psFree(file->fpa); file->fpa = NULL; - if (data->masksName) { - pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.MASK", j); - psFree(file->fpa); - file->fpa = NULL; - } pmFPAfileActivate(data->config->files, false, NULL); } @@ -373,8 +343,4 @@ 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"); - } } @@ -389,6 +355,4 @@ ro1->image = image1; ro2->image = image2; - ro1->mask = mask1; - ro2->mask = mask2; ro1->data_exists = cell1->data_exists = cell1->parent->data_exists = true; @@ -409,17 +373,11 @@ 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; @@ -437,20 +395,8 @@ 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 } } @@ -458,19 +404,7 @@ 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); -#if WCS_DEBUG - 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 - ); -#endif modify_wcs1 = 0; } @@ -482,12 +416,16 @@ pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1 -/* // 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); - psMetadataAddF64(Fcell1->concepts,PS_LIST_TAIL,"EXPTIME", PS_META_REPLACE, "", stats->robustMedian); + + // I am baffled that this is the way to get the exposure time updated correctly. + psMetadataItem *item1; + item1 = psMetadataLookup(Fcell1->concepts,"CELL.EXPOSURE"); + item1->data.F32 = stats->robustMedian; + pmReadout *Fro1 = pmReadoutAlloc(Fcell1); Fro1->image = image1; - Fro1->mask = mask1; Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; @@ -495,31 +433,4 @@ 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); @@ -527,5 +438,5 @@ fits2->fpa->hdu->header = psMetadataAlloc(); psMetadataCopy(fits2->fpa->hdu->header,projhdu->header); - + if (modify_wcs2) { pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header); @@ -561,9 +472,10 @@ psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); - psMetadataAddF64(Fcell2->concepts,PS_LIST_TAIL,"EXPTIME", PS_META_REPLACE, "", stats->robustMedian); + + psMetadataItem *item2 = psMetadataLookup(Fcell2->concepts,"CELL.EXPOSURE"); + item2->data.F32 = stats->robustMedian; pmReadout *Fro2 = pmReadoutAlloc(Fcell2); Fro2->image = image2; - Fro2->mask = mask2; Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true; @@ -575,23 +487,4 @@ psFree(view); -#if 0 - { - psString filename = NULL; // Filename for image - psStringAppend(&filename, "skycell_%d_1.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, image1, 0, NULL); - psFitsClose(fits); - } - - { - psString filename = NULL; // Filename for image - psStringAppend(&filename, "skycell_%d_2.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, image2, 0, NULL); - psFitsClose(fits); - } -#endif filesIterateUp(data->config);