Index: trunk/psphot/src/psModulesUtils.c
===================================================================
--- trunk/psphot/src/psModulesUtils.c	(revision 6379)
+++ trunk/psphot/src/psModulesUtils.c	(revision 6481)
@@ -144,8 +144,10 @@
     int xIs, xJs, yIs, yJs;
     int xIe, yIe;
-    float flux;
+    float flux, wt;
 
     psImage *Pi = Mi->pixels;
     psImage *Pj = Mj->pixels;
+
+    psImage *Wi = Mi->weight;
 
     psImage *Ti = Mi->mask;
@@ -166,4 +168,6 @@
     yIe = Ye - Pi->row0;
 
+    // note that this is addressing the same image pixels,
+    // though only if both are source not model images
     flux = 0;
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
@@ -171,5 +175,54 @@
 	    if (Ti->data.U8[yi][xi]) continue;
 	    if (Tj->data.U8[yj][xj]) continue;
-	    flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
+	    wt = Wi->data.F32[yi][xi];
+	    if (wt > 0) {
+		flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+	    } 
+	}
+    }
+    return (flux);
+}
+
+float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    float flux, wt;
+
+    psImage *Pi = Mi->pixels;
+    psImage *Pj = Mj->pixels;
+
+    psImage *Wi = Mi->weight;
+
+    psImage *Ti = Mi->mask;
+    psImage *Tj = Mj->mask;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+    
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+    
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that this is addressing the same image pixels,
+    // though only if both are source not model images
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+	for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+	    if (Ti->data.U8[yi][xi]) continue;
+	    if (Tj->data.U8[yj][xj]) continue;
+	    wt = Wi->data.F32[yi][xi];
+	    if (wt > 0) {
+		flux += 1.0 / wt;
+	    } 
 	}
     }
@@ -250,2 +303,11 @@
     return NULL;
 }
+
+bool psRegionIsNaN (psRegion region) {
+
+    if (!isfinite(region.x0)) return true;
+    if (!isfinite(region.x1)) return true;
+    if (!isfinite(region.y0)) return true;
+    if (!isfinite(region.y1)) return true;
+    return false;
+}