IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

new version of the burntool programme

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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   }
Note: See TracChangeset for help on using the changeset viewer.