Changeset 27840 for branches/simtest_nebulous_branches/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (modified) (9 diffs)

branches/simtest_nebulous_branches/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c

@@ -25 +25 @@
    err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell);
 
+//   fprintf(stderr, "Got through burn_check\n");
+
 /* Expand the masks; burns asymmetrically, stars symmetrically */
 //   err = grow_mask(nx, ny, NX, mbuf, YMASK_GROW, YMASK_GROW, cell->nstar, cell->star);
@@ -30 +32 @@
    err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, 
                    MASK_SAT_HALO, cell->nburn, cell->burn);
+//   fprintf(stderr, "Got through grow mask for burns\n");
    err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, 
                    MASK_STAR_HALO, cell->nstar, cell->star);
+//   fprintf(stderr, "Got through grow mask for stars\n");
 
 /* Fix up all the burns */
@@ -37 +41 @@
       if(!cell->burn[k].burned) continue;
 /* Fit the trail */
+//      fprintf(stderr, "Fitting trail %d\n", k);
       fit_trail(nx, ny, NX, imbuf, mbuf, cell->burn+k, 1, 
              cell->sky+cell->bias, cell->rms, BURN_PWR);
 /* Subtract the fit from the image */
+//      fprintf(stderr, "Subtracting trail %d\n", k);
       if(!cell->burn[k].fiterr) sub_fit(nx, ny, NX, buf, cell->burn+k, 1);
    }
@@ -103 +109 @@
 /* Restore all the burns */
    for(k=0; k<cell->npersist; k++) {
-      if(!cell->persist[k].fiterr) 
+      if(!cell->persist[k].fiterr && 
+         (cell->persist[k].func == BURN_PWR || 
+          cell->persist[k].func == BURN_EXP) ) {
          sub_fit(nx, ny, NX, buf, &(cell->persist[k]), -1);
+      }
    }
    return(0);
@@ -118 +127 @@
 /* Restore all the burns */
    for(k=0; k<cell->npersist; k++) {
-      if(!cell->persist[k].fiterr) 
+      if(!cell->persist[k].fiterr && 
+         (cell->persist[k].func == BURN_PWR || 
+          cell->persist[k].func == BURN_EXP) ) {
          sub_fit(nx, ny, NX, buf, &(cell->persist[k]), +1);
+      }
    }
    return(0);
@@ -218 +230 @@
 /****************************************************************/
 /* burn_test() tests whether a box is burned or not */
+/* It's trying to balance the flux of above-below to see whether 
+ * there's a net burn above.  It has to cope with the case that the box is
+ * so low there is no above in which case it compares center-side.
+ * It's somewhat rudimentary as indicated by the FIXME's, and the original
+ * version had some real logical problems.
+ */
 STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms,
                        MTYPE *mask, OBJBOX *box)
 {
+   int i, j, nburn, nref, nmed;
+   int y0, y1, ymid, xburn, xref, dx=0, dy=0;
+
+/* FIXME: needs a test for flat-toppedness as well as trail... */
+/* FIXME: needs a test for pure, massive saturation */
+
+/* Center line */
+   ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+
+/* Irrelevant, too near the top to detect a burn anyway */
+   if(box->ey > ny-FIT_EDGE-3) {
+      box->diff = -1;
+      box->burned = -1;
+/* But looks like a bad one which will leave behind persistence */
+      if(box->ey-box->sy+1 > 4) box->burned = 1;
+      return(0);
+
+/* Too near the bottom, do a side-side test */
+   } else if(box->sy < FIT_EDGE+3) {
+      dy = 0;
+      y0 = ymid + MAX(box->ey-ymid, FIT_EDGE);
+      y1 = MIN(y0+100, ny-1);
+
+      if(2*box->ex - box->sx + 1 < nx-1) {      /* Work right? */
+         dx = box->ex - box->sx + 1;
+
+      } else if(2*box->sx - box->ex - 1 > 0) {  /* Work left? */
+         dx = -(box->ex - box->sx + 1);
+
+      } else {  /* No room!  Oh no! */
+         box->diff = -1;
+         box->burned = -1;
+/* But looks like a bad one which will leave behind persistence */
+         if(box->ey-box->sy+1 > 4) box->burned = 1;
+         return(0);
+      }
+
+/* Do a top-bottom test for excess flux */
+   } else {
+      dx = 0;
+      dy = MAX(ymid-box->sy, box->ey-ymid);
+      if(dy < FIT_EDGE) dy = FIT_EDGE;
+      y0 = ymid + dy;
+      y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1);
+      if(y1-y0 > 100) y1 = y0 + 100;
+   }
+
+/* Do the test; get the median of the averages across x */
+   nmed = 0;
+   for(j=y0; j<y1; j++) {
+      xburn = xref = 0;
+      nburn = nref = 0;
+      for(i=box->sx; i<=box->ex; i++) {
+         if(mask[i+j*NX] == MASK_NONE) {
+            xburn += data[i+j*NX];
+            nburn++;
+         }
+         if(dx != 0) {  /* Side to side */
+            if(mask[i+dx+j*NX] == MASK_NONE) {
+               xref += data[i+dx+j*NX];
+               nref++;
+            }
+         } else {       /* top-bottom */
+            if(mask[i+(2*ymid-j)*NX] == MASK_NONE) {
+               xref += data[i+(2*ymid-j)*NX];
+               nref++;
+            }
+         }
+      }
+      if(nburn > 0) xburn /= nburn;
+      if(nref > 0) xref /= nref;
+      if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref;
+   }
+   box->diff = int_median(nmed, median_buf);
+   box->burned = (box->diff > 0.3*rms);
+   return(0);
+}
+
+#ifdef ORIG_BURN_TEST
+// This has a lot of problems -- chief among which is what happens with
+// a huge kite-shaped bad region.  It's fundamentally trying to balance
+// the flux of above-below to see whether there's a net burn above.  But
+// it has to cope with the case that the box is so low there is no above
+// in which case it compares center-side.
+/****************************************************************/
+/* burn_test() tests whether a box is burned or not */
+STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms,
+                       MTYPE *mask, OBJBOX *box)
+{
    int i, j, n, nrow, nmed;
    int y0, y1, ymid, xsum, dx=0;
@@ -229 +336 @@
 /* Center line */
    ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+
 /* Starting point offset */
    y0 = MAX(ymid-box->sy, box->ey-ymid);
@@ -259 +367 @@
       nrow = 0;
       for(i=box->sx; i<=box->ex; i++) {
+// IS THIS THE PROBLEM?
          xsum += data[i+(ymid+y0+j)*NX]*(mask[i+(ymid+y0+j)*NX] == MASK_NONE);
          if(y1 < ny) {
@@ -279 +388 @@
    return(0);
 }
+#endif