IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42392 for trunk


Ignore:
Timestamp:
Feb 14, 2023, 12:48:54 PM (3 years ago)
Author:
tdeboer
Message:

new version of the burntool programme

Location:
trunk/extsrc/gpcsw/gpcsrc/fits/burntool
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c

    r26766 r42392  
    2323/* Find all the burned patches and bright stars */
    2424   err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum);
     25
    2526   err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell);
    2627
    … …  
    3435//   fprintf(stderr, "Got through grow mask for burns\n");
    3536   err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW,
    36                    MASK_STAR_HALO, cell->nstar, cell->star);
     37         MASK_STAR_HALO, cell->nstar, cell->star);
    3738//   fprintf(stderr, "Got through grow mask for stars\n");
    3839
    … …  
    7071   printf("  #    cx  cy  max     sx  sy    ex  ey  midy  diff  S/B/E  slope     zero\n");
    7172   for(k=0; k<cell->nburn; k++) {
    72       ymid = (cell->burn[k].y0m + cell->burn[k].y1m +
    73               cell->burn[k].y0p + cell->burn[k].y1p + 2) / 4;
     73      /*ymid = (cell->burn[k].y0m + cell->burn[k].y1m +
     74              cell->burn[k].y0p + cell->burn[k].y1p + 2) / 4;*/
     75      ymid = (cell->burn[k].sy + cell->burn[k].ey + 1) / 2;
    7476      i = cell->burn[k].nfit / 2;
    7577      printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %5d %2d %1d %1d %6.3f %8d\n",
    … …  
    8890   printf("  #    cx  cy  max     sx  sy    ex  ey  midy  S/P\n");
    8991   for(k=0; k<cell->nstar; k++) {
    90       ymid = (cell->star[k].y0m + cell->star[k].y1m +
    91               cell->star[k].y0p + cell->star[k].y1p + 2) / 4;
     92      /*ymid = (cell->star[k].y0m + cell->star[k].y1m +
     93              cell->star[k].y0p + cell->star[k].y1p + 2) / 4;*/
     94      ymid = (cell->star[k].sy + cell->star[k].ey + 1) / 2;
    9295
    9396      printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %1d\n",
    … …  
    142145                       MTYPE *mask, CELL *cell)
    143146{
    144    int err, k;
     147   int err,k;
    145148
    146149/* Given these big stars, identify which ones really have a burn */
    … …  
    160163             k, cell->burn[k].sx, cell->burn[k].sy,
    161164             cell->burn[k].ex, cell->burn[k].ey,
    162              cell->burn[k].y0m, cell->burn[k].y1m, cell->burn[k].y0p, cell->burn[k].y0p,
     165             cell->burn[k].y0m, cell->burn[k].y1m, cell->burn[k].y0p, cell->burn[k].y1p,
    163166             cell->burn[k].diff, cell->burn[k].burned);
    164167   }
    … …  
    242245   int y0, y1, ymid, xburn, xref, dx=0, dy=0;
    243246
     247   /*introduce a maximum to the burn contrast check, to avoid a very large number in crowded regions*/
     248   int burn_contrast = MIN(0.3*rms,5);
     249   /*introduce a maximum size for the burn contrast check box, instead of a random large 100 pixels */
     250   /*for small stars you will lose signal after 10 pixels, so use the fitted box size as a good guide*/
     251   int test_box = MAX(box->ey-box->sy,10);
     252
     253   if(VERBOSE & VERB_BOXGROW) {
     254        printf("burn contrast test using box of max size: %8d\n", test_box);
     255   }
     256
    244257/* FIXME: needs a test for flat-toppedness as well as trail... */
    245258/* FIXME: needs a test for pure, massive saturation */
    246259
    247260/* Center line */
    248    ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
     261   //ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
     262   ymid = (box->sy+box->ey)/2;
    249263
    250264/* Irrelevant, too near the top to detect a burn anyway */
    … …  
    283297      y0 = ymid + dy;
    284298      y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1);
    285       if(y1-y0 > 100) y1 = y0 + 100;
     299      if(y1-y0 > test_box) y1 = y0 + test_box;
    286300   }
    287301
    … …  
    292306      nburn = nref = 0;
    293307      for(i=box->sx; i<=box->ex; i++) {
    294          if(mask[i+j*NX] == MASK_NONE) {
    295             xburn += data[i+j*NX];
    296             nburn++;
    297          }
    298          if(dx != 0) {  /* Side to side */
    299             if(mask[i+dx+j*NX] == MASK_NONE) {
    300                xref += data[i+dx+j*NX];
    301                nref++;
    302             }
    303          } else {       /* top-bottom */
    304             if(mask[i+(2*ymid-j)*NX] == MASK_NONE) {
    305                xref += data[i+(2*ymid-j)*NX];
    306                nref++;
    307             }
    308          }
     308             if(mask[i+j*NX] == MASK_NONE) {
     309               xburn += data[i+j*NX];
     310               nburn++;
     311             }
     312             if(dx != 0) {      /* Side to side */
     313               if(mask[i+dx+j*NX] == MASK_NONE) {
     314                xref += data[i+dx+j*NX];
     315                nref++;
     316               }
     317             } else {   /* top-bottom */
     318               if(mask[i+(2*ymid-j)*NX] == MASK_NONE) {
     319                xref += data[i+(2*ymid-j)*NX];
     320                nref++;
     321               }
     322             }
    309323      }
    310324      if(nburn > 0) xburn /= nburn;
    311325      if(nref > 0) xref /= nref;
    312326      if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref;
     327      if(VERBOSE & VERB_BOXGROW) {
     328        printf("contrast: %8d %8d %8d %8d %8d\n", j, xburn, (2*ymid-j), xref, xburn - xref);
     329      }
     330
    313331   }
    314332   box->diff = int_median(nmed, median_buf);
    315    box->burned = (box->diff > 0.3*rms);
     333   box->burned = (box->diff > burn_contrast);
    316334   return(0);
    317335}
    … …  
    335353
    336354/* Center line */
    337    ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
     355   //ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
     356   ymid = (box->sy + box->ey + 1) / 2;
    338357
    339358/* Starting point offset */

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

    r26766 r42392  
    8989   PSF_THRESH   =  5000;        /* Threshold for a star to be a PSF */
    9090
    91    BMASK_GROW =  5.0;           /* Growth of burned boxes in size */
    92    RMASK_GROW =  3.0;           /* Growth of burn/star in diameter */
     91   BMASK_GROW =  3.0;           /* Growth of burned boxes in pixels */
     92   RMASK_GROW =  2.0;           /* Growth of burn/star in diameter */
    9393
    9494   MIN_PSF_SIZE  =    3;        /* Min box size for stamp selection */
    … …  
    9797   MAX_PSF_PER_CELL = 20;       /* Max number of PSF stars accepted per cell */
    9898
    99    NEGLIGIBLE_TRAIL = 0.5;      /* Don't sweat less than this * sigma */
     99   NEGLIGIBLE_TRAIL = 0.4;      /* Don't sweat less than this * sigma */
    100100   EXPIRE_TRAIL_TIME = 2000;    /* Expire a persist after this [sec] */
    101101
    … …  
    564564
    565565/* Does this cell look kosher? */
    566          if(OTA[cell].rms*OTA[cell].rms >
    567             OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
    568             if(VERBOSE > 0) {
    569                fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
    570                fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n",
    571                        OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
    572             }
    573             continue;
    574          }
     566//       if(OTA[cell].rms*OTA[cell].rms >
     567//          OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
     568//          if(VERBOSE > 0) {
     569//             fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
     570//             fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n",
     571//                     OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
     572//          }
     573//          continue;
     574//       }
    575575
    576576         if(apply) {
    … …  
    601601            if(VERBOSE & VERB_NORM) burn_blab(OTA+cell);
    602602
    603 /* Fix up the streaks (don't bother if table only) */
    604 // Why the heck do I get to skip this if table only?  Makes no sense!
    605 //          if(!tableonly) {
    606                persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf,
     603            persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf,
    607604                           OTA+cell);
    608605/* Tell us about it? */
    609                if(VERBOSE & VERB_NORM) persist_blab(OTA+cell);
    610 //          }
     606            if(VERBOSE & VERB_NORM) persist_blab(OTA+cell);
     607
    611608         }
    612609
    … …  
    667664{
    668665/* Make some space for medians */
    669    if(2*SKY_MARG*ny > nmedian_buf) {
     666/*   if(2*SKY_MARG*ny > nmedian_buf) {
    670667      if(median_buf != NULL) free(median_buf);
    671668      if( (median_buf = (int *)calloc(2*SKY_MARG*ny, sizeof(int))) == NULL) {
    … …  
    675672
    676673      nmedian_buf = 2*SKY_MARG*ny;
     674   }
     675*/
     676
     677/* Make some space for cell copy */
     678   if(nx*ny > nmedian_buf) {
     679      if(median_buf != NULL) free(median_buf);
     680      if( (median_buf = (int *)calloc(nx*ny, sizeof(DTYPE))) == NULL) {
     681       fprintf(stderr, "\rerror: failed to alloc cell copy\n");
     682       exit(-668);
     683      }
     684      nmedian_buf = nx*ny;
    677685   }
    678686
    … …  
    728736/* Get sky stats */
    729737   for(k=n=nsat=0; k<nx*ny; k+=((617*nx)/1000)) {
     738//   for(k=n=nsat=0; k<nx*ny; k+=((100*nx)/1000)) {
     739//   for(k=n=nsat=0; k<nx*ny; k++) {
    730740      i = k % nx;
    731741      j = k / nx;

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h

    r26766 r42392  
    7070#define STAR_RADIUS  4          /* Radius over which a star ctr must be max */
    7171#define SKY_MARG 3              /* Horiz offset from burn area for sky */
    72 #define FIT_EDGE 10             /* How far beyond saturation to start fit? */
     72#define FIT_EDGE 5              /* How far beyond saturation to start fit? */
    7373#define Y_SCALE 0.01            /* Scale factor for y in fits */
    7474
    … …  
    102102
    103103/* Fit parameters */
    104 #define FIT_MIN_SLOPE  -5.0     /* minimum slope which is a credible fit */
     104#define FIT_MIN_SLOPE -10.0     /* minimum slope which is a credible fit */
    105105#define FIT_MAX_SLOPE   0.0     /* maximum slope which is a credible fit */
    106106
    … …  
    144144      int sxfit;        /* starting column for fits */
    145145      int exfit;        /* ending column for fits */
     146      int eyfit;  /* y-coord ending column for fits */
    146147      int nfit;         /* how many columns were corrected? */
    147148      IMTYPE *stamp;    /* postage stamp of this object */

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

    r26766 r42392  
    1818STATIC int grow_box(int nx, int ny, int NX, DTYPE *data, int thr, OBJBOX *box)
    1919{
    20    int i, j, is, ie, left, right;
     20   int i, j,i2,j2, is, ie,jlo,jhi,ilo,ihi, left, right;
    2121
    2222   left = right = 1;
    … …  
    2727
    2828   do {
    29 /* Push the right wall until it's clear */
     29      /* Push the right wall until it's clear */
     30      jlo = box->sy;
     31      jhi = box->ey;
    3032      if(right) {
    31          for(i=box->ex+1; i<nx; i++) {
    32             for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
    33             if(j > box->ey) break;
    34             if(VERBOSE & VERB_BOXGROW) {
    35                printf("right: %8d %8d\n", i, j);
    36             }
    37          }
    38          box->ex = i - 1;
    39          right = 0;
    40       }
    41 
    42 /* Push the left wall until it's clear */
     33             for(i=box->ex+1; i<nx; i++) {
     34               for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
     35          if(j > box->ey) break;
     36          for(j2=box->ey; j2>=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
     37          if(j2 < box->sy) break;
     38
     39          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
     40          /*If not, there is a discontinuity, and we might have overlapping sources*/
     41          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
     42
     43               if(VERBOSE & VERB_BOXGROW) {
     44                 printf("right: %8d %8d %8d\n", i, j,data[i+j*NX]);
     45               }
     46          jlo = j;
     47          jhi = j2;
     48             }
     49             box->ex = i - 1;
     50             right = 0;
     51      }
     52
     53      /* Push the left wall until it's clear */
     54      jlo = box->sy;
     55      jhi = box->ey;
    4356      if(left) {
    44          for(i=box->sx-1; i>=0; i--) {
    45             for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
    46             if(j > box->ey) break;
    47             if(VERBOSE & VERB_BOXGROW) {
    48                printf("left: %8d %8d\n", i, j);
    49             }
    50          }
    51          box->sx = i + 1;
    52          left = 0;
    53       }
    54 
    55 /* Push the top wall until it's clear, catching corners! */
     57             for(i=box->sx-1; i>=0; i--) {
     58               for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
     59               if(j > box->ey) break;
     60          for(j2=box->ey; j2<=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
     61          if(j2 < box->sy) break;
     62
     63          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
     64          /*If not, there is a discontinuity, and we might have overlapping sources*/
     65          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
     66
     67               if(VERBOSE & VERB_BOXGROW) {
     68                 printf("left: %8d %8d %8d\n", i, j,data[i+j*NX]);
     69               }
     70          jlo = j;
     71          jhi = j2;
     72             }
     73             box->sx = i + 1;
     74             left = 0;
     75      }
     76
     77      /* Push the top wall until it's clear, catching corners! */
    5678      is = MAX(0, box->sx-1);
    5779      ie = MIN(nx-1, box->ex+1);
     80      ilo = is;
     81      ihi = ie;
    5882      for(j=box->ey+1; j<ny; j++) {
    59          for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
    60          if(i == is && box->sx>0) left = 1;     /* Need to redo left */
    61          if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
    62          if(i > ie) break;                      /* All clear */
    63          if(VERBOSE & VERB_BOXGROW) {
    64             printf("top: %8d %8d\n", i, j);
    65          }
     83             for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
     84        if(i > ie) break;        /* All clear */
     85        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
     86        if(i2 < is) break;        /* All clear */
     87
     88             if(i == is && box->sx>0) left = 1; /* First point hits thr, so need to redo left */
     89             if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /*Last point hits thr, so redo right */
     90
     91        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
     92        /*If not, there is a discontinuity, and we might have overlapping sources*/
     93        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) {break;}
     94
     95             if(VERBOSE & VERB_BOXGROW) {
     96               printf("top: %8d %8d %8d\n", i, j,data[i+j*NX]);
     97             }
     98        ilo = i;
     99        ihi = i2;
    66100      }
    67101      box->ey = j - 1;
    68102
    69 /* Push the bottom wall until it's clear, catching corners! */
     103      /* Push the bottom wall until it's clear, catching corners! */
    70104      is = MAX(0,box->sx-1);
    71105      ie = MIN(nx-1,box->ex+1);
     106      ilo = is;
     107      ihi = ie;
    72108      for(j=box->sy-1; j>=0; j--) {
    73          for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
    74          if(i == is && box->sx>0) left = 1;     /* Need to redo left */
    75          if(data[is+1+j*NX] >= thr) box->y0m = j;
    76          if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
    77          if(data[ie-1+j*NX] >= thr) box->y1m = j;
    78          if(i > ie) break;                      /* All clear */
    79          if(VERBOSE & VERB_BOXGROW) {
    80             printf("bottom: %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m);
    81          }
     109             for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
     110        if(i > ie) break;        /* All clear */
     111        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
     112        if(i2 < is) break;        /* All clear */
     113
     114             if(i == is && box->sx>0) left = 1; /* Need to redo left */
     115        if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
     116
     117             if(data[is+1+j*NX] >= thr) box->y0m = j;
     118             if(data[ie-1+j*NX] >= thr) box->y1m = j;
     119
     120        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
     121        /*If not, there is a discontinuity, and we might have overlapping sources*/
     122        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) break;
     123
     124             if(VERBOSE & VERB_BOXGROW) {
     125               printf("bottom: %8d %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m,data[i+j*NX]);
     126             }
     127        ilo = i;
     128        ihi = i2;
    82129      }
    83130      box->sy = j + 1;
    … …  
    126173                     double rfac, int maskval, int nbox, OBJBOX *box)
    127174{
    128    int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1;
     175   int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1,boxsize;
    129176
    130177   for(k=0; k<nbox; k++) {
    … …  
    137184
    138185      for(j=j0; j<=j1; j++) {
    139          for(i=i0+1; i<i1; i++) {
    140             if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
    141             if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
    142          }
     186             for(i=i0+1; i<i1; i++) {
     187              if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
     188              if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
     189             }
    143190      }
    144191
    145192/* Push out the center to a diameter rfac*xsize */
    146193      x0 = (box[k].sx + box[k].ex + 1) / 2;
    147       y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
    148       dr = NINT((box[k].ex - box[k].sx + 1) * 0.5 * rfac);
     194      //y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
     195      y0 = (box[k].sy + box[k].ey + 1) / 2;
     196      boxsize = MAX((box[k].ex - box[k].sx + 1),(box[k].ey - box[k].sy + 1));
     197      dr = NINT(boxsize * 0.5 * rfac);
    149198      for(j=MAX(0,y0-dr); j<=MIN(ny-1, y0+dr); j++) {
    150          for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
    151             r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
    152             if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
     199             for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
     200               r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
     201               if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
    153202               mask[i+j*NX] = maskval;
    154          }
     203             }
    155204      }
    156205   }

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c

    r26766 r42392  
    2323
    2424/* Merge all overlapping streaks */
    25    persist_merge(cell);
     25  persist_merge(cell);
     26
    2627
    2728/* Fix up all the persistence streaks */
    2829   for(k=0; k<cell->npersist; k++) {
    29 
    30 /* Is this just a blasted area being carried for IPP? */
    31       if((cell->persist)[k].func == BURN_BLASTED) {
    32          if(cell->time - (cell->persist)[k].time > EXPIRE_TRAIL_TIME) {
     30     /* Is this just a blasted area being carried for IPP? */
     31     if((cell->persist)[k].func == BURN_BLASTED) {
     32           if(cell->time - (cell->persist)[k].time > EXPIRE_TRAIL_TIME) {
    3333            (cell->persist)[k].fiterr = FIT_EXPIRED;
    34          }
    35          continue;
    36       }
     34           }
     35           continue;
     36     }
    3737
    3838/* This had a significant positive slope, what do we do now? */
    39       if( (cell->persist)[k].func == BURN_POSSLOPE) {
     39     if( (cell->persist)[k].func == BURN_POSSLOPE) {
    4040/* try again, so let it slide through... */
    41       }
     41     }
    4242
    4343/* Fit the trail */
    44       err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0,
     44     err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0,
    4545             cell->sky+cell->bias, cell->rms, BURN_EXP);
    4646
    4747/* Subtract out the fit */
    48       if(!cell->persist[k].fiterr)
    49          err = sub_fit(nx, ny, NX, buf, cell->persist+k, 1);
     48     if(!cell->persist[k].fiterr)
     49           err = sub_fit(nx, ny, NX, buf, cell->persist+k, 1);
    5050   }
    5151
    … …  
    6565   for(k=0; k<cell->npersist; k++) {
    6666      i = (cell->persist[k].ex - cell->persist[k].sx + 1) / 2;
    67       ymid = (cell->persist[k].y0m + cell->persist[k].y1m +
    68               cell->persist[k].y0p + cell->persist[k].y1p + 2) / 4;
    69       printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %8.5f %4d\n",
     67      /*ymid = (cell->persist[k].y0m + cell->persist[k].y1m +
     68              cell->persist[k].y0p + cell->persist[k].y1p + 2) / 4;*/
     69        ymid = (cell->persist[k].sy + cell->persist[k].ey + 1) / 2;
     70      printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %8.5f %8.5f\n",
    7071             k, cell->persist[k].cx, cell->persist[k].cy,
    7172             cell->persist[k].max,
    … …  
    7374             cell->persist[k].ex, cell->persist[k].ey, ymid,
    7475             cell->persist[k].func,
    75              cell->persist[k].slope, NINT((cell->persist[k]).zero[i]));
     76             cell->persist[k].slope, (cell->persist[k]).zero[i]);
    7677   }
    7778   return(0);

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c

    r37678 r42392  
    3838   while(fgets(line, 1024, fp) != NULL) {
    3939      if(line[0] == '#') continue;
    40       sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d %d",
     40      sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d %d %d",
    4141             &boxbuf[nbox].cell, &boxbuf[nbox].time,
    4242             &boxbuf[nbox].cx, &boxbuf[nbox].cy,
    … …  
    5151             &boxbuf[nbox].slope, &boxbuf[nbox].nfit,
    5252             &boxbuf[nbox].sxfit, &boxbuf[nbox].exfit,
    53              &boxbuf[nbox].fiterr);
     53             &boxbuf[nbox].fiterr, &boxbuf[nbox].eyfit);
    5454      if(boxbuf[nbox].nfit > 0) {
    55          boxbuf[nbox].zero = (double *)calloc(boxbuf[nbox].nfit, sizeof(double));
    56          boxbuf[nbox].xfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
    57          boxbuf[nbox].yfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
    58          if(boxbuf[nbox].zero == NULL ||
    59             boxbuf[nbox].xfit == NULL ||
    60             boxbuf[nbox].yfit == NULL) {
    61             fprintf(stderr, "\rerror: failed to alloc boxbuf\n");
    62             exit(-673);
    63          }
     55             boxbuf[nbox].zero = (double *)calloc(boxbuf[nbox].nfit, sizeof(double));
     56             boxbuf[nbox].xfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
     57             boxbuf[nbox].yfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
     58             if(boxbuf[nbox].zero == NULL ||
     59              boxbuf[nbox].xfit == NULL ||
     60              boxbuf[nbox].yfit == NULL) {
     61             fprintf(stderr, "\rerror: failed to alloc boxbuf\n");
     62             exit(-673);
     63           }
    6464
    6565         for(i=0; i<boxbuf[nbox].nfit; i++) {
    … …  
    123123{
    124124   int i, j, k, n, nbox, xs, xe;
    125    int jp, kp;
     125   int jp, kp,xdiff;
    126126   int lapmax;
    127127   OBJBOX *box;
    … …  
    137137   for(i=0; i<myMAXSIZE; i++) boxid[i] = -1;
    138138   for(k=0; k<nbox; k++) {
    139       if(box[k].exfit >= myMAXSIZE-1) continue;
    140       for(i=box[k].sxfit; i<=box[k].exfit; i++) xusage[i] += 1;
     139      if(box[k].ex >= myMAXSIZE-1) continue;
     140      for(i=box[k].sx; i<=box[k].ex; i++) xusage[i] += 1;
    141141   }
    142142
    143143/* Identify clusters */
    144144   for(xs=0; xs<myMAXSIZE-1; xs++) {
    145       if(xusage[xs] == 0) continue;
     145    if(xusage[xs] == 0) continue;
    146146      for(xe=xs, lapmax=0; xe<myMAXSIZE-1; xe++) {
    147          if(xusage[xe+1] == 0) break;
    148          lapmax = MAX(lapmax, xusage[xe]);
     147             if(xusage[xe+1] == 0) break;
     148             lapmax = MAX(lapmax, xusage[xe]);
    149149      }
    150150      if(lapmax == 1) { /* No overlap?  No problem. */
    151          xs = xe + 1;   /* Hop to next gap */
    152          continue;
    153       }
    154 
    155 /* Which boxes overlap this cluster? */
     151            xs = xe + 1;        /* Hop to next gap */
     152            continue;
     153      }
     154
     155      /* Which boxes overlap this cluster? */
    156156      for(k=n=0; k<nbox; k++) {
    157          if(box[k].sxfit > xe || box[k].exfit < xs) continue;
    158          boxid[n] = k;
    159          yctr[n] = (box[k].y0m+box[k].y1m+box[k].y0p+box[k].y0p+2) / 4;
    160          n++;
    161       }
    162 
    163 /* Case 1: different y start => sever */
     157        if(box[k].sx > xe || box[k].ex < xs) continue;
     158             boxid[n] = k;
     159             //yctr[n] = (box[k].y0m+box[k].y1m+box[k].y0p+box[k].y1p+2) / 4;
     160             yctr[n] = (box[k].sy+box[k].ey+1) / 2;
     161             n++;
     162      }
     163
    164164      for(kp=0; kp<n; kp++) {
    165          k = boxid[kp];
    166          zk = 0.0;
    167          if(box[k].nfit > 0) zk = box[k].zero[box[k].nfit/2];
    168          for(jp=kp+1; jp<n; jp++) {
    169             j = boxid[jp];
    170             zj = 0.0;
    171             if(box[j].nfit > 0) zj = box[j].zero[box[j].nfit/2];
    172             if(ABS(yctr[jp]-yctr[kp]) > DIFFERENT_STREAK) {
    173 /* Trim back the feebler streak */
    174                if(zk > zj) {
    175                   if(box[j].sxfit >= box[k].sxfit)
    176                      box[j].sxfit = MAX(box[j].sxfit, box[k].exfit+1);
    177                   if(box[j].exfit <= box[k].exfit)
    178                      box[j].exfit = MIN(box[j].exfit, box[k].sxfit-1);
    179                } else {
    180                   if(box[k].sxfit >= box[j].sxfit)
    181                      box[k].sxfit = MAX(box[k].sxfit, box[j].exfit+1);
    182                   if(box[k].exfit <= box[j].exfit)
    183                      box[k].exfit = MIN(box[k].exfit, box[j].sxfit-1);
    184                }
    185             }
    186          }
    187       }
    188 
    189 /* Case 2: pretty much the same y start => union */
    190       for(kp=0; kp<n; kp++) {
    191          k = boxid[kp];
    192          for(jp=kp+1; jp<n; jp++) {
    193             j = boxid[jp];
    194             /* CZW: Since I added the initialization statement above, */
    195             /*      any box with boxid == -1 hasn't been set, and */
    196             /*      therefore needs to be dropped, I think. */
    197             if ((k < 0)||(j < 0)) { continue; }
    198             if(ABS(yctr[jp]-yctr[kp]) <= DIFFERENT_STREAK) {
    199 /* Merge k into j, trash k from further consideration */
     165            k = boxid[kp];
     166            zk = 0.0;
     167            if(box[k].nfit > 0) zk = box[k].zero[box[k].nfit/2];
     168            for(jp=kp+1; jp<n; jp++) {
     169             j = boxid[jp];
     170             zj = 0.0;
     171             if(box[j].nfit > 0) zj = box[j].zero[box[j].nfit/2];
     172             if(box[j].sx > box[k].sx) xdiff = box[j].sx - box[k].ex;
     173             if(box[j].sx <= box[k].sx) xdiff = box[k].sx - box[j].ex;
     174
     175             /* CZW: Since I added the initialization statement above, */
     176             /*      any box with boxid == -1 hasn't been set, and */
     177             /*      therefore needs to be dropped, I think. */
     178             if ((k < 0)||(j < 0)) { continue; }
     179
     180        /*Tdb20220222: In the case of overlapping x-coords but not overlapping y-coords*/
     181        /*I fail to see the need to change the x-range of any streaks. That just spells disaster*/
     182        // Case 1: different y start => sever
     183             if(ABS(yctr[jp]-yctr[kp]) > DIFFERENT_STREAK) {
     184         // Trim back the feebler streak
     185              if(zk > zj) {
     186                    /*
     187                    if(box[j].sxfit >= box[k].sxfit)
     188               box[j].sxfit = MAX(box[j].sxfit, box[k].exfit+1);
     189                    if(box[j].exfit <= box[k].exfit)
     190                    box[j].exfit = MIN(box[j].exfit, box[k].sxfit-1);
     191                    */
     192              } else {
     193               /*       
     194          if(box[k].sxfit >= box[j].sxfit)
     195                    box[k].sxfit = MAX(box[k].sxfit, box[j].exfit+1);
     196                    if(box[k].exfit <= box[j].exfit)
     197                    box[k].exfit = MIN(box[k].exfit, box[j].sxfit-1);
     198                    */
     199              }
     200             }
     201
     202             /* Case 2: pretty much the same y start => union */
     203             if((ABS(yctr[jp]-yctr[kp]) <= DIFFERENT_STREAK) && (xdiff < 0) ) {
     204          /* Merge k into j, trash k from further consideration */
    200205               box[j].time = MAX(box[j].time, box[k].time);
    201206               box[j].sx = MIN(box[j].sx, box[k].sx);
    … …  
    213218               box[j].sxfit = MIN(box[j].sxfit, box[k].sxfit);
    214219               box[j].exfit = MAX(box[j].exfit, box[k].exfit);
    215 
    216                box[k].exfit = box[k].sxfit - 1;
     220               box[j].eyfit = MIN(box[j].eyfit, box[k].eyfit);
     221
     222               box[k].exfit = -999;
    217223               yctr[kp] = -2 * DIFFERENT_STREAK;
    218224            }
    … …  
    226232/* Excise the boxes which have been eradicated (sxfit > exfit) */
    227233   for(k=0; k<nbox; k++) {
    228       if(box[k].sxfit > box[k].exfit) {
    229          for(j=k; j<nbox-1; j++) memcpy(box+j, box+j+1, sizeof(OBJBOX));
    230          nbox--;
    231          k--;
     234      if(box[k].exfit < 0) {
     235            for(j=k; j<nbox-1; j++) memcpy(box+j, box+j+1, sizeof(OBJBOX));
     236            nbox--;
     237            k--;
    232238      }
    233239   }
    … …  
    256262   fprintf(fp,  "# Cell: %d  sky= %d   rms= %d   bias= %d\n",
    257263           cell[0].cell, cell[0].sky, cell[0].rms, cell[0].bias);
    258    fprintf(fp, "#Cell time    cx  cy  max   y0   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p F up    slope nfit sxf exf\n");
     264   fprintf(fp, "#Cell time    cx  cy  max   y0   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p F up    slope nfit sxf exf fiterr eyf\n");
    259265
    260266   for(j=0; j<MAXCELL; j++) {
    … …  
    276282            if(cell[j].persist[k].nfit <= 0) continue;
    277283         }
    278          fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d\n",
     284
     285         fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d %3d\n",
    279286                 j, cell[j].persist[k].time,
    280287                 cell[j].persist[k].cx, cell[j].persist[k].cy,
    … …  
    289296                 cell[j].persist[k].slope, cell[j].persist[k].nfit,
    290297                 cell[j].persist[k].sxfit, cell[j].persist[k].exfit,
    291                  cell[j].persist[k].fiterr);
     298                 cell[j].persist[k].fiterr, cell[j].persist[k].eyfit);
    292299         for(i=0; i<cell[j].persist[k].nfit; i++) {
    293300            fprintf(fp, "%3d %3d %8.4f\n", cell[j].persist[k].xfit[i],
    … …  
    298305/* Second: new burns */
    299306      for(k=0; k<cell[j].nburn; k++) {
    300          if(!cell[j].burn[k].burned) continue;
    301          if(PERSIST_RETAIN) {
    302 /* Keep fits which have a dubious slope */
    303             if(cell[j].burn[k].fiterr != FIT_SLOPE_ERROR) {
     307            if(!cell[j].burn[k].burned) continue;
     308            if(PERSIST_RETAIN) {
     309        /* Keep fits which have a dubious slope */
     310             if(cell[j].burn[k].fiterr != FIT_SLOPE_ERROR) {
    304311               if(cell[j].burn[k].fiterr &&
    305                   cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
     312                    cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
    306313               if(cell[j].burn[k].nfit <= 0) continue;
     314             }
     315            } else {
     316             if(cell[j].burn[k].fiterr &&
     317               cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
     318             if(cell[j].burn[k].nfit <= 0) continue;
    307319            }
    308          } else {
    309             if(cell[j].burn[k].fiterr &&
    310                cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
    311             if(cell[j].burn[k].nfit <= 0) continue;
    312          }
    313320
    314321         i = (cell[j].burn[k].ex - cell[j].burn[k].sx + 1) / 2;
    315          fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d\n",
     322         fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d %3d\n",
    316323                 j, cell[j].burn[k].time,
    317324                 cell[j].burn[k].cx, cell[j].burn[k].cy,
    … …  
    326333                 cell[j].burn[k].slope, cell[j].burn[k].nfit,
    327334                 cell[j].burn[k].sxfit, cell[j].burn[k].exfit,
    328                  cell[j].burn[k].fiterr);
     335                 cell[j].burn[k].fiterr, cell[j].burn[k].eyfit);
    329336         for(i=0; i<cell[j].burn[k].nfit; i++) {
    330337            fprintf(fp, "%3d %3d %8.4f\n", cell[j].burn[k].xfit[i],

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c

    r26766 r42392  
    3434   starcut = STAR_THRESH/2 + cell->sky + cell->bias;
    3535
     36   /*try a new value for thresh_lo which is 1 sigmas above the background, but with a safety limit of 1000 counts*/
     37        trailthresh = cell->sky + (10*cell->rms);
     38        if(trailthresh <= 1000) trailthresh = 1000;
     39        trailthresh = trailthresh + cell->bias;
     40
    3641/* Look at all the pixels which pass the burn threshold */
    3742   for(j=0; j<ny; j++) {
    3843      xon = -1;
    3944      for(i=0; i<nx; i++) {
    40 /* Big enough?  Initialize a box and push it outwards to the trail level. */
    41          if(mask[i+j*NX] > MASK_NONE) continue;
    42          if(xon < 0) {
    43             if(data[i+j*NX] > burnthresh) {
    44                xon = i;
    45                thresh_hi = burnthresh;
    46                thresh_lo = trailthresh;
    47                if(VERBOSE & VERB_DETECT) {
    48                   printf("Starting a burn at %d %d %d\n", i, j, data[i+j*NX]);
    49                }
    50             } else if(data[i+j*NX] > starthresh) {
    51 /* Local maximum? */
    52                if(!local_max(i, j, STAR_RADIUS, nx, ny, NX, data)) continue;
    53                if(data[i+j*NX] < veto[i+j*NX]) continue;
    54                xon = i;
    55 //             thresh_hi = starthresh;
    56 //             thresh_lo = starcut;
    57                thresh_hi = data[i+j*NX];
    58                thresh_lo = STAR_FRAC*(data[i+j*NX]-cell->sky-cell->bias) +
    59                   cell->sky + cell->bias;
    60                if(VERBOSE & VERB_DETECT) {
    61                   printf("Starting a star at %d %d %d\n", i, j, data[i+j*NX]);
    62                }
    63             }
    64          }
    65 
    66          if(xon >= 0 && (data[i+j*NX] < thresh_hi || i == nx-1)) {
    67 
    68 //       if(xon < 0 && data[i+j*NX] > burnthresh && mask[i+j*NX] == 0) xon = i;
    69 //       if(xon > 0 && (data[i+j*NX] < burnthresh || i == nx-1)) {
    70 
    71             xoff = i<nx-1 ? i-1 : nx-1;
    72             if(nbox >= MAXBURN) {
    73                fprintf(stderr, "error: too many burn boxes\n");
    74                return(-1);
    75             }
    76             boxbuf[nbox].sx = xon;
    77             boxbuf[nbox].ex = xoff;
    78             boxbuf[nbox].sy = j;
    79             for(k=j; k<ny && data[(xon+xoff)/2+k*NX] > thresh_hi; k++);
    80             boxbuf[nbox].ey = k-1;
    81             grow_box(nx, ny, NX, data, thresh_lo, boxbuf+nbox);
    82 /* Fill in max and center info */
    83             boxbuf[nbox].max = 0;
    84             for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
    85                for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
    86                   if(data[k+l*NX] > boxbuf[nbox].max) {
    87                      boxbuf[nbox].cx = k;
    88                      boxbuf[nbox].cy = l;
    89                      boxbuf[nbox].max = data[k+l*NX];
    90                   }
    91                }
    92             }
     45         /* Big enough?  Initialize a box and push it outto trail level. */
     46              if(mask[i+j*NX] > MASK_NONE) continue;
     47              if(xon < 0) {
     48                if(data[i+j*NX] > burnthresh) {
     49                  xon = i;
     50                  thresh_hi = burnthresh;
     51                  thresh_lo = trailthresh;
     52                  if(VERBOSE & VERB_DETECT) {
     53                         printf("Starting a burn at %d %d %d %d\n", i, j, data[i+j*NX],thresh_lo);
     54                  }
     55                } else if(data[i+j*NX] > starthresh) {
     56             /* Local maximum? */
     57                  if(!local_max(i, j, STAR_RADIUS, nx, ny, NX, data)) continue;
     58                  if(data[i+j*NX] < veto[i+j*NX]) continue;
     59                  xon = i;
     60//                thresh_hi = starthresh;
     61//                thresh_lo = starcut;
     62                  thresh_hi = data[i+j*NX];
     63                  thresh_lo = STAR_FRAC*(data[i+j*NX]-cell->sky-cell->bias) +
     64                       cell->sky + cell->bias;
     65                  if(VERBOSE & VERB_DETECT) {
     66                         printf("Starting a star at %d %d %d %d\n", i, j, data[i+j*NX],thresh_lo);
     67                  }
     68                }
     69              }
     70
     71              if(xon >= 0 && (data[i+j*NX] < thresh_hi || i == nx-1)) {
     72                xoff = i<nx-1 ? i-1 : nx-1;
     73                if(nbox >= MAXBURN) {
     74                   fprintf(stderr, "error: too many burn boxes\n");
     75                   return(-1);
     76                }
     77                boxbuf[nbox].sx = xon;
     78                boxbuf[nbox].ex = xoff;
     79                boxbuf[nbox].sy = j;
     80                for(k=j; k<ny && data[(xon+xoff)/2+k*NX] > thresh_hi; k++);
     81                boxbuf[nbox].ey = k-1;
     82                grow_box(nx, ny, NX, data, thresh_lo, boxbuf+nbox);
     83           /* Fill in max and center info */
     84                boxbuf[nbox].max = 0;
     85                for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
     86                 for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
     87                       if(data[k+l*NX] > boxbuf[nbox].max) {
     88                        boxbuf[nbox].cx = k;
     89                        boxbuf[nbox].cy = l;
     90                        boxbuf[nbox].max = data[k+l*NX];
     91                       }
     92                 }
     93                }
     94
    9395/* A box which triggered on starthresh cannot envelop burnthresh */
    94             if(thresh_hi < burnthresh && boxbuf[nbox].max > burnthresh) {
    95                if(VERBOSE & VERB_DETECT) {
    96                   printf("Ditching box %d %d %d %d  max = %d > burnthresh\n",
    97                          boxbuf[nbox].sx, boxbuf[nbox].sy,
    98                          boxbuf[nbox].ex, boxbuf[nbox].ey,
    99                          boxbuf[nbox].max);
    100                   fflush(stdout);
    101                }
     96                if(thresh_hi < burnthresh && boxbuf[nbox].max > burnthresh) {
     97                  if(VERBOSE & VERB_DETECT) {
     98                        printf("Ditching box %d %d %d %d  max = %d > burnthresh\n",
     99                             boxbuf[nbox].sx, boxbuf[nbox].sy,
     100                             boxbuf[nbox].ex, boxbuf[nbox].ey,
     101                             boxbuf[nbox].max);
     102                        fflush(stdout);
     103                  }
    102104/* Mark the veto mask so as not to trigger on this one again */
    103                for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
    104                   for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
    105                      veto[k+l*NX] = boxbuf[nbox].max;
    106                   }
    107                }
    108                xon = -1;
    109                continue;
    110             }
     105                  for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
     106                         for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
     107                          veto[k+l*NX] = boxbuf[nbox].max;
     108                         }
     109                  }
     110                  xon = -1;
     111                  continue;
     112                }
    111113
    112114/* Mask the pixels so as not to catch them again! */
    113             for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
    114                for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
    115                   mask[k+l*NX] = MASK_CTR;
    116                }
    117             }
    118             boxbuf[nbox].sat = boxbuf[nbox].max >= burnthresh;
    119             boxbuf[nbox].max -= cell->sky + cell->bias;
     115                for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
     116                 for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
     117                       mask[k+l*NX] = MASK_CTR;
     118                 }
     119                }
     120                boxbuf[nbox].sat = boxbuf[nbox].max >= burnthresh;
     121                boxbuf[nbox].max -= cell->sky + cell->bias;
    120122/* Burned only if it triggered on burnthresh, not if it enveloped it */
    121             boxbuf[nbox].burned = thresh_hi >= burnthresh;
     123                boxbuf[nbox].burned = thresh_hi >= burnthresh;
    122124/* But take a hard look at stars with really bright centers... */
    123             if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) {
    124                burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox);
    125             }
    126             boxbuf[nbox].time = cell->time;
    127             boxbuf[nbox].cell = cellnum;
    128             boxbuf[nbox].stamp = NULL;
    129             boxbuf[nbox].func = FUNC_NONE;
    130             boxbuf[nbox].sxfit = boxbuf[nbox].sx;
    131             boxbuf[nbox].exfit = boxbuf[nbox].ex;
    132             boxbuf[nbox].nfit = 0;
    133             boxbuf[nbox].zero = NULL;
    134             boxbuf[nbox].xfit = boxbuf[nbox].yfit = NULL;
    135             if(VERBOSE & VERB_DETECT) {
    136                printf("New box at %d %d %d %d",
     125                if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) {
     126                 burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox);
     127                }
     128                boxbuf[nbox].time = cell->time;
     129                boxbuf[nbox].cell = cellnum;
     130                boxbuf[nbox].stamp = NULL;
     131                boxbuf[nbox].func = FUNC_NONE;
     132 
     133                boxbuf[nbox].sxfit = boxbuf[nbox].sx;
     134                boxbuf[nbox].exfit = boxbuf[nbox].ex;
     135                boxbuf[nbox].nfit = 0;
     136                boxbuf[nbox].zero = NULL;
     137                boxbuf[nbox].xfit = boxbuf[nbox].yfit = NULL;
     138                if(VERBOSE & VERB_DETECT) {
     139                 printf("New box at %d %d %d %d",
    137140                      boxbuf[nbox].sx, boxbuf[nbox].sy,boxbuf[nbox].ex, boxbuf[nbox].ey);
    138                printf("   %d %d %d %d",
     141                 printf("   %d %d %d %d",
    139142                      boxbuf[nbox].x0m, boxbuf[nbox].x0p,boxbuf[nbox].x1m, boxbuf[nbox].x1p);
    140                printf("   %d %d %d %d\n",
     143                 printf("   %d %d %d %d\n",
    141144                      boxbuf[nbox].y0m, boxbuf[nbox].y0p,boxbuf[nbox].y1m, boxbuf[nbox].y1p);
    142                fflush(stdout);
    143             }
    144             xon = -1;
    145             nbox++;
    146          }
     145                 fflush(stdout);
     146                }
     147                xon = -1;
     148                nbox++;
     149              }
    147150      }
    148151   }

  • trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c

    r26766 r42392  
    1414#include "burnparams.h"
    1515
    16 STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE];
     16STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE];
    1717
    1818/****************************************************************/
    … …  
    8686#define MAXWGT 2.0      /* Maximum allowed weight */
    8787#define MINWGT 0.1      /* Minimum allowed weight */
    88 #define MAXDEV 0.5      /* Maximum ln deviation above first pass fit */
     88#define MAXDEV 1.       /* Maximum ln deviation above first pass fit */
    8989#define MAXRMS 4.0      /* Minimum factor of RMS to be chopped */
    9090
    … …  
    9494                    int up, int bckgnd, int rms, int fitfunc)
    9595{
    96    int i, j, k, err, nfit;
     96   int i, j, k, err, nfit,firstfitpix;
    9797   int xs, xe, y0, y1, y2, dy, yfit;
    98    double slope, zero, zsum, wsum, trial;
     98   double slope, zero, zsum, wsum, trial,zfun;
    9999   char fooname[60];
    100100   FILE *fp=stdout;
    … …  
    105105   }
    106106
    107    xs = box->sxfit;
    108    xe = box->exfit;
    109 
     107   /*xs = box->sxfit; */
     108   /*xe = box->exfit; */
     109   /*TdB20220222: Do not use s/exfit here, which will fail in the case of positive slopes*/
     110   xs = box->sx;
     111   xe = box->ex;
     112   
    110113   if(up) {
    111       y0 = (box->y0m + box->y1m + box->y0p + box->y0p + 2) / 4;
    112 /* Maybe a bit better?  Maybe not... */
    113 //      y0 = box->ey;
    114       y1 = box->ey + FIT_EDGE;  /* Start of fit, budged away from the burn */
    115       y2 = ny - 1;              /* End of fit */
     114      //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
     115      y0 = (box->sy + box->ey + 1) / 2;
     116      y1 = box->ey + FIT_EDGE;   /* Start of fit, budged away from the burn */
     117      y2 = ny - 1;      /* End of fit */
    116118      dy = +1;
    117       yfit = box->sy;           /* Start of validity, TBD later by col */
     119      yfit = box->sy;      /* Start of validity, TBD later by col */
     120      firstfitpix = MAX((box->ey - box->sy)*2,20) ;  /* Number of pixels to include in the initial fit*/
    118121   } else {
    119       y0 = (box->y0m + box->y1m + box->y0p + box->y0p + 2) / 4;
    120       y1 = 0;                   /* Start of fit */
    121       y2 = box->sy - FIT_EDGE;  /* End of fit, budged away from the burn */
    122 //      y2 = box->ey - FIT_EDGE;        /* End of fit, budged away from the burn */
    123 //      y2 = y0;        /* End of fit, budged away from the burn */
     122      //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
     123      y0 = (box->sy + box->ey + 1) / 2;
     124      y1 = 0;        /* Start of fit */
     125      y2 = y0 - FIT_EDGE;  /* End of fit, budged away from the burn */
    124126      dy = -1;
    125       yfit = box->ey;           /* Start of validity, TBD later by col */
     127      yfit = box->ey;      /* Start of validity, TBD later by col */
     128      firstfitpix = 100;  /* Number of pixels to include in the initial fit*/
     129
    126130   }
    127131
    … …  
    147151   box->up = up;
    148152   box->y0 = y0;
     153   box->eyfit = up?y2:y1;
    149154   for(i=0; i<xe-xs+1; i++) {
    150155      box->zero[i] = 0.0;
    … …  
    157162      box->slope = -1.0;
    158163      box->fiterr = FIT_TOP_ERROR;
    159 /* No fit here because the burn is to top, but the persist may be fitable */
     164      /* No fit here because the burn is to top, but the persist may be fitable */
    160165      if(box->sy > FIT_EDGE) {
    161          box->nfit = xe - xs + 1;
    162 /* However, if burn extends to bottom, persist will not be able to fit */
     166            box->nfit = xe - xs + 1;
     167       /* However, if burn extends to bottom, persist will not be able to fit */
    163168      } else {
    164          box->nfit = 1;         /* Just enough to keep it alive */
    165          box->func = BURN_BLASTED;
     169            box->nfit = 1;              /* Just enough to keep it alive */
     170            box->func = BURN_BLASTED;
    166171      }
    167172      return(0);
    … …  
    174179
    175180/* Accumulate points to fit */
     181/* use wbuf0 to select only the first FIRSTFITPIX pixels for first fit */
    176182   nfit = y2 - y1 + 1;
    177183   for(j=y1; j<=y2; j++) {
    … …  
    179185      zbuf[j]=0.0;
    180186      for(i=xs, k=0; i<=xe; i++) {
    181 //       printf("%4d %4d %7d %7d\n", i, j, data[i+j*NX], mask[i+j*NX]);
    182          if(data[i+j*NX] != NODATA && mask[i+j*NX] == MASK_NONE) {
    183 //       if(data[i+j*NX] != NODATA) {
    184             zbuf[j] += data[i+j*NX];
    185             k++;
    186          }
    187       }
    188 //      printf("%4d %4d %7d\n", k, j, zbuf[j]);
     187             if(data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO) ) {
     188              zbuf[j] += data[i+j*NX];
     189              k++;
     190             }
     191        if(mask[i+j*NX] != MASK_NONE) {
     192        }
     193      }
    189194      zbuf[j] = zbuf[j] / MAX(1,k) - bckgnd;
    190 /* Linearize for fit... */
     195
     196/*    Linearize for fit... */
    191197      if(fitfunc == BURN_PWR) ybuf[j] = log(ybuf[j]);
    192198      wbuf[j] = 0;
     199      wbuf0[j] = 0;
    193200      if(zbuf[j] > 1) {
    194          wbuf[j] = MIN(MAXWGT, (zbuf[j] / rms));  /* Quadratic not good */
    195          wbuf[j] = MAX(MINWGT, wbuf[j]);
    196          zbuf[j] = log(zbuf[j]);
     201             wbuf[j] = MIN(MAXWGT, (zbuf[j] / rms));  /* Quadratic not good */
     202             wbuf[j] = MAX(MINWGT, wbuf[j]);
     203             zbuf[j] = log(zbuf[j]);
     204        wbuf0[j] = wbuf[j];
     205        if(j < (y2-firstfitpix) && !up) wbuf0[j] = 0;
     206        if(j > (y1+firstfitpix) && up) wbuf0[j] = 0;
     207        //printf("fit: %4d %.3f %.3f %.3f %.3f\n", j, ybuf[j], zbuf[j], wbuf[j], wbuf0[j]);
    197208      }
    198209   }
    199210/* First pass fit */
    200    err = wlinearfit(nfit, ybuf+y1, zbuf+y1, wbuf+y1, &slope, &zero);
     211   err = wlinearfit(nfit, ybuf+y1, zbuf+y1, wbuf0+y1, &slope, &zero);
    201212   if(VERBOSE & VERB_FIT) {
    202213      printf("nfit= %d  err= %d  slope= %.3f  zero= %.3f\n",
    … …  
    210221/* Trim away any really big positive deviations */
    211222   for(j=y1; j<=y2; j++) {
    212       if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV &&
    213          zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
     223//      if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV &&
     224//       zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
     225      if(abs(zbuf[j]-(ybuf[j]*slope+zero)) > MAXDEV) wbuf[j] = 0;
    214226      if(VERBOSE & VERB_FITPROF) {
    215227         fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j);
    … …  
    239251//      if(trial > 2*rms) box->nfit = -box->nfit;
    240252      if(trial > 2*rms) box->func = BURN_POSSLOPE;
     253      box->sxfit = xs;
     254      box->exfit = xs - 1;
    241255      return(-1);
    242256   }
    … …  
    244258/* Reconstruct the fit for comparison with the data col by col */
    245259   for(j=y1; j<=y2; j++) zbuf[j] = exp(ybuf[j]*slope);
    246 /* And tag on the extras to determine the end of the fit */
     260   /* And tag on the extras to determine the end of the fit */
    247261   for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
    248262      if(fitfunc == BURN_EXP) {
    249 //       zbuf[j] = exp(Y_SCALE*MAX(0,dy*(j-y0))*slope);
    250          zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
     263             zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
    251264      } else {
    252          if(dy*(j-y0) <= 0) zbuf[j] = exp(log(Y_SCALE*1.0)*slope);
    253          else zbuf[j] = exp(log(Y_SCALE*dy*(j-y0))*slope);
    254       }
    255    }
    256 
    257 /* Save the fit information */
     265             if(dy*(j-y0) <= 0) zbuf[j] = exp(log(Y_SCALE*1.0)*slope);
     266             else zbuf[j] = exp(log(Y_SCALE*dy*(j-y0))*slope);
     267      }
     268   }
     269
     270/* find where the fitted function becomes less than half the rms and cut the trail there*/
     271   for(j=y1; j<=y2; j++) {
     272      zfun = exp((ybuf[j]*slope)+zero);
     273      if(zfun > NEGLIGIBLE_TRAIL*rms) break;
     274   }   
     275   //For burns, just let it go to the top
     276   if(j < y2 && !up) box->eyfit = j;
     277
     278/* Determine good x-ranges and Save the fit information */
    258279   box->slope = slope;
    259280   box->nfit = 0;
    … …  
    262283      zsum = wsum = 0.0;
    263284      for(j=y1; j<=y2; j++) {
    264          if(wbuf[j] > 0 && data[i+j*NX] != NODATA &&
    265             mask[i+j*NX] == MASK_NONE) {
    266 //       if(wbuf[j] > 0 && data[i+j*NX] != NODATA) {
    267 //          wsum += wbuf[j];
    268 //          zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j];
    269             wsum += zbuf[j];
    270             zsum += data[i+j*NX] - bckgnd;
    271          }
     285             if(wbuf[j] > 0 && data[i+j*NX] != NODATA &&
     286             (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) {
     287//             if(wbuf[j] > 0 && data[i+j*NX] != NODATA) {
     288//             wsum += wbuf[j];
     289//             zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j];
     290               wsum += zbuf[j];
     291               zsum += data[i+j*NX] - bckgnd;
     292             }
    272293      }
    273294      if(zsum < 0 || wsum <= 0) {
    274          zsum = 0.0;
     295             zsum = 0.0;
    275296      } else {
    276          zsum = zsum / wsum;
    277       }
    278 /* FIXME: what's a really good criterion for negligible fit? */
    279 /* 100 pixels up fit is zsum or ~zsum/e */
    280 //      if(zsum > NEGLIGIBLE_TRAIL*rms) {
     297             zsum = zsum / wsum;
     298      }
     299      /* FIXME: what's a really good criterion for negligible fit? */
     300      /* 100 pixels up fit is zsum or ~zsum/e */
     301//    if(zsum > NEGLIGIBLE_TRAIL*rms) {
     302      /* 100113: but also evaluate at the end of the fit for stubby trails */
    281303      if(zsum > NEGLIGIBLE_TRAIL*rms ||
    282 /* 100113: but also evaluate at the end of the fit for stubby trails */
    283          zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
    284 
    285 /* Ascertain the starting point of where the fit is good */
     304           zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
     305        /* Ascertain the starting point of where the fit is good */
    286306//      box->yfit[box->nfit] = yfit;    /* Vanilla, misses center */
    287307//      box->yfit[box->nfit] = y0;      /* Center, no adjust */
    288          for(j=yfit; dy*j<=dy*(up?y2:y1); j+=dy) {
    289             trial = zsum * zbuf[j];
    290 /* First time data - fit is closer to sky than is data */
    291             if(data[i+j*NX] > bckgnd + 0.5*trial) break;
    292          }
    293          box->yfit[box->nfit] = j;
    294          box->zero[box->nfit] = zsum;
    295          box->xfit[box->nfit] = i;
    296          box->nfit += 1;
     308             for(j=yfit; dy*j<=dy*(up?y2:y1); j+=dy) {
     309               trial = zsum * zbuf[j];
     310/*        First time data - fit is closer to sky than is data */
     311               if(data[i+j*NX] > bckgnd + 0.5*trial) break;
     312             }
     313             box->yfit[box->nfit] = j;
     314             box->zero[box->nfit] = zsum;
     315             box->xfit[box->nfit] = i;
     316             box->nfit += 1;
    297317      }
    298318   }
    … …  
    306326      box->fiterr = FIT_ALL_GONE;
    307327   }
     328
    308329//   if(box->nfit <= 0) fprintf(stderr, "WHOA, got nfit = 0\n");
    309330   return(0);
Note: See TracChangeset for help on using the changeset viewer.