Index: trunk/ppSkycell/src/ppSkycellLoop.c
===================================================================
--- trunk/ppSkycell/src/ppSkycellLoop.c	(revision 34800)
+++ trunk/ppSkycell/src/ppSkycellLoop.c	(revision 35852)
@@ -10,4 +10,7 @@
 
 #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
@@ -22,4 +25,5 @@
 }
 
+
 static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System
                                int numCols, int numRows // Size of image
@@ -34,13 +38,14 @@
     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; \
+    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);
@@ -49,7 +54,8 @@
     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
@@ -150,9 +156,17 @@
         pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System
         if (!wcs) {
+	  if (data->wcsrefName) {
+	    pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", i);
+	    wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
+	  }
+	  if (!wcs) {	  
             psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", i);
             return false;
-        }
-
-        psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image
+	  }
+        }
+
+	psRegion *region = psRegionAlloc(-1.0 * wcs->crpix1,numCols - wcs->crpix1,
+					 -1.0 * wcs->crpix2,numRows - wcs->crpix2);
+	imageRegions->data[i] = region;
         psTrace("ppSkycell", 5, "Image region %d is: [%.0f:%.0f,%.0f:%.0f]\n",
                 i, region->x0, region->x1, region->y0, region->y1);
@@ -162,5 +176,10 @@
             if (wcs->crval1 == crval1->data.F64[j] && wcs->crval2 == crval2->data.F64[j] &&
                 wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) {
-                regionMinMax(projRegions->data[j], region);
+	      //                regionMinMax(projRegions->data[j], region);
+	      psRegion *proj = projRegions->data[j];
+	      proj->x0 = PS_MIN(region->x0,proj->x0);
+	      proj->x1 = PS_MAX(region->x1,proj->x1);
+	      proj->y0 = PS_MIN(region->y0,proj->y0);
+	      proj->y1 = PS_MAX(region->y1,proj->y1);
                 target->data.S32[i] = j;
                 found = true;
@@ -188,7 +207,11 @@
     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 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
@@ -209,7 +232,15 @@
             psImageInit(mask2, 0xFF);
         }
+	// HDU containing the WCS we plan on using
 	pmHDU *projhdu = NULL;
+
+	// Do we need to modify the WCS?  This flips to zero after we've done it once.
 	int modify_wcs1 = 1;
 	int modify_wcs2 = 1;
+
+	// Because we may have holes, we need to ensure that we set the CRPIX in the binned image correction.
+	// Find the minimum/maximum, so we know where the zero is.
+	float maxCRPIX1   = -99e99;
+	float maxCRPIX2   = -99e99;
 	// Loop over inputs to this projection.
         for (int j = 0; j < data->numInputs; j++) {
@@ -229,25 +260,36 @@
             }
 
+	    // Read the HDU/WCS information from the first entry, and use that as the reference.
             pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
+	    pmAstromWCS *wcs = pmAstromWCSfromHeader(file->fpa->hdu->header); // World Coordinate System
+	    if (!wcs) {
+	      if (data->wcsrefName) {
+		pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", j);
+		wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
+		if (!projhdu) {
+		  projhdu = wcsref->fpa->hdu;
+		}
+	      }
+	      if (!wcs) {	  
+		psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", j);
+		return false;
+	      }
+	    }
 	    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	  
-	    
+	    }
+
+	    // However, we need to check to see if we've found the maximum CRPIX
+	    if (wcs->crpix1 > maxCRPIX1) {
+	      maxCRPIX1 = wcs->crpix1;
+	    }
+	    if (wcs->crpix2 < maxCRPIX2) {
+	      maxCRPIX2 = wcs->crpix2;
+	    }
+
+	    // Actuall do the binning.
             pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
             psFree(view);
 
-	    //	    data->maskVal = 0xffff;
             pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
             if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
@@ -262,9 +304,22 @@
             }
 
+#if (0)
             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->y1 + projRegion->y1) / (float)data->bin1;
+	    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;
+	    psTrace("ppSkycell",5,"Offsets: %d %d : %d %d",
+		    xOffset1,yOffset1,xOffset2,yOffset2);
+
+	    // 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
@@ -276,4 +331,5 @@
             }
 
+	    // Cleanup on input loop.
             psFree(bin1RO);
             psFree(bin2RO);
@@ -345,6 +401,17 @@
 	    WCS->cdelt1 *= cd1f;
 	    WCS->cdelt2 *= cd2f;
-	    WCS->crpix1 = WCS->crpix1 / cd1f;
-	    WCS->crpix2 = WCS->crpix2 / cd2f;
+	    // Fudge the CRPIX incase we have missing corners
+	    if (maxCRPIX1 > WCS->crpix1) {
+	      WCS->crpix1 = maxCRPIX1 / cd1f;
+	    }
+	    else {
+	      WCS->crpix1 = WCS->crpix1 / cd1f;
+	    }
+	    if (maxCRPIX2 > WCS->crpix2) {
+	      WCS->crpix2 = maxCRPIX2 / cd2f;
+	    }
+	    else {
+	      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,
@@ -376,6 +443,6 @@
 	    }
 	    pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
+#if WCS_DEBUG
 	    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,
@@ -448,11 +515,11 @@
 	    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++) {
+	    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++) {
+	    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);
 	      }