Changeset 42392 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c

Timestamp:

Feb 14, 2023, 12:48:54 PM (3 years ago)

Author:

tdeboer

Message:

new version of the burntool programme

File:

• trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c (modified) (4 diffs)

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c

@@ -18 +18 @@
 STATIC int grow_box(int nx, int ny, int NX, DTYPE *data, int thr, OBJBOX *box)
 {
-   int i, j, is, ie, left, right;
+   int i, j,i2,j2, is, ie,jlo,jhi,ilo,ihi, left, right;
 
    left = right = 1;
@@ -27 +27 @@
 
    do {
-/* Push the right wall until it's clear */
+      /* Push the right wall until it's clear */
+      jlo = box->sy;
+      jhi = box->ey;
       if(right) {
-         for(i=box->ex+1; i<nx; i++) {
-            for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
-            if(j > box->ey) break;
-            if(VERBOSE & VERB_BOXGROW) {
-               printf("right: %8d %8d\n", i, j);
-            }
-         }
-         box->ex = i - 1;
-         right = 0;
-      }
-
-/* Push the left wall until it's clear */
+             for(i=box->ex+1; i<nx; i++) {
+               for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
+          if(j > box->ey) break;
+          for(j2=box->ey; j2>=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
+          if(j2 < box->sy) break;
+
+          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+          /*If not, there is a discontinuity, and we might have overlapping sources*/
+          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
+
+               if(VERBOSE & VERB_BOXGROW) {
+                 printf("right: %8d %8d %8d\n", i, j,data[i+j*NX]);
+               }
+          jlo = j;
+          jhi = j2;
+             }
+             box->ex = i - 1;
+             right = 0;
+      }
+
+      /* Push the left wall until it's clear */
+      jlo = box->sy;
+      jhi = box->ey;
       if(left) {
-         for(i=box->sx-1; i>=0; i--) {
-            for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
-            if(j > box->ey) break;
-            if(VERBOSE & VERB_BOXGROW) {
-               printf("left: %8d %8d\n", i, j);
-            }
-         }
-         box->sx = i + 1;
-         left = 0;
-      }
-
-/* Push the top wall until it's clear, catching corners! */
+             for(i=box->sx-1; i>=0; i--) {
+               for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
+               if(j > box->ey) break;
+          for(j2=box->ey; j2<=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
+          if(j2 < box->sy) break;
+
+          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+          /*If not, there is a discontinuity, and we might have overlapping sources*/
+          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
+
+               if(VERBOSE & VERB_BOXGROW) {
+                 printf("left: %8d %8d %8d\n", i, j,data[i+j*NX]);
+               }
+          jlo = j;
+          jhi = j2;
+             }
+             box->sx = i + 1;
+             left = 0;
+      }
+
+      /* Push the top wall until it's clear, catching corners! */
       is = MAX(0, box->sx-1);
       ie = MIN(nx-1, box->ex+1);
+      ilo = is;
+      ihi = ie;
       for(j=box->ey+1; j<ny; j++) {
-         for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
-         if(i == is && box->sx>0) left = 1;     /* Need to redo left */
-         if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
-         if(i > ie) break;                      /* All clear */
-         if(VERBOSE & VERB_BOXGROW) {
-            printf("top: %8d %8d\n", i, j);
-         }
+             for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
+        if(i > ie) break;        /* All clear */
+        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
+        if(i2 < is) break;        /* All clear */
+
+             if(i == is && box->sx>0) left = 1; /* First point hits thr, so need to redo left */
+             if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /*Last point hits thr, so redo right */
+
+        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+        /*If not, there is a discontinuity, and we might have overlapping sources*/
+        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) {break;}
+
+             if(VERBOSE & VERB_BOXGROW) {
+               printf("top: %8d %8d %8d\n", i, j,data[i+j*NX]);
+             }
+        ilo = i;
+        ihi = i2;
       }
       box->ey = j - 1;
 
-/* Push the bottom wall until it's clear, catching corners! */
+      /* Push the bottom wall until it's clear, catching corners! */
       is = MAX(0,box->sx-1);
       ie = MIN(nx-1,box->ex+1);
+      ilo = is;
+      ihi = ie;
       for(j=box->sy-1; j>=0; j--) {
-         for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
-         if(i == is && box->sx>0) left = 1;     /* Need to redo left */
-         if(data[is+1+j*NX] >= thr) box->y0m = j;
-         if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
-         if(data[ie-1+j*NX] >= thr) box->y1m = j;
-         if(i > ie) break;                      /* All clear */
-         if(VERBOSE & VERB_BOXGROW) {
-            printf("bottom: %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m);
-         }
+             for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
+        if(i > ie) break;        /* All clear */
+        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
+        if(i2 < is) break;        /* All clear */
+
+             if(i == is && box->sx>0) left = 1; /* Need to redo left */
+        if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
+
+             if(data[is+1+j*NX] >= thr) box->y0m = j;
+             if(data[ie-1+j*NX] >= thr) box->y1m = j;
+
+        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+        /*If not, there is a discontinuity, and we might have overlapping sources*/
+        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) break;
+
+             if(VERBOSE & VERB_BOXGROW) {
+               printf("bottom: %8d %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m,data[i+j*NX]);
+             }
+        ilo = i;
+        ihi = i2;
       }
       box->sy = j + 1;
@@ -126 +173 @@
                      double rfac, int maskval, int nbox, OBJBOX *box)
 {
-   int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1;
+   int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1,boxsize;
 
    for(k=0; k<nbox; k++) {
@@ -137 +184 @@
 
       for(j=j0; j<=j1; j++) {
-         for(i=i0+1; i<i1; i++) {
-            if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
-            if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
-         }
+             for(i=i0+1; i<i1; i++) {
+              if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
+              if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
+             }
       }
 
 /* Push out the center to a diameter rfac*xsize */
       x0 = (box[k].sx + box[k].ex + 1) / 2;
-      y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
-      dr = NINT((box[k].ex - box[k].sx + 1) * 0.5 * rfac);
+      //y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
+      y0 = (box[k].sy + box[k].ey + 1) / 2;
+      boxsize = MAX((box[k].ex - box[k].sx + 1),(box[k].ey - box[k].sy + 1));
+      dr = NINT(boxsize * 0.5 * rfac);
       for(j=MAX(0,y0-dr); j<=MIN(ny-1, y0+dr); j++) {
-         for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
-            r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
-            if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
+             for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
+               r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
+               if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
                mask[i+j*NX] = maskval;
-         }
+             }
       }
    }