IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42753


Ignore:
Timestamp:
Dec 6, 2024, 1:28:47 PM (20 months ago)
Author:
eugene
Message:

updates from trunk: burntool improvements

Location:
branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c

    r42392 r42753  
    1717/* burn_fix(): Find and fix all the burned patches in a cell */
    1818STATIC int burn_fix(int nx, int ny, int NX, int NY, IMTYPE *buf,
    19                     CELL *cell, int cellnum)
     19                    CELL *cell, int cellnum, char *camera)
    2020{
    2121   int k, err;
    2222
    2323/* Find all the burned patches and bright stars */
    24    err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum);
    25 
    26    err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell);
     24   err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum, camera);
     25
     26   err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell, camera);
    2727
    2828//   fprintf(stderr, "Got through burn_check\n");
    2929
    3030/* Expand the masks; burns asymmetrically, stars symmetrically */
    31 //   err = grow_mask(nx, ny, NX, mbuf, YMASK_GROW, YMASK_GROW, cell->nstar, cell->star);
    32 //   err = grow_mask(nx, ny, NX, mbuf, XMASK_GROW, YMASK_GROW, cell->nburn, cell->burn);
    3331   err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW,
    3432                   MASK_SAT_HALO, cell->nburn, cell->burn);
    … …  
    4442//      fprintf(stderr, "Fitting trail %d\n", k);
    4543      fit_trail(nx, ny, NX, imbuf, mbuf, cell->burn+k, 1,
    46              cell->sky+cell->bias, cell->rms, BURN_PWR);
     44             cell->sky+cell->bias, cell->rms, BURN_PWR, camera);
    4745/* Subtract the fit from the image */
    4846//      fprintf(stderr, "Subtracting trail %d\n", k);
    … …  
    5452      for(k=0; k<NX*ny; k++) buf[k] = mbuf[k] - BZERO;
    5553   }
     54   
    5655   return(err);
    5756}
    … …  
    106105/****************************************************************/
    107106/* burn_restore(): Restore the fitted trails back to the image */
    108 STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell)
     107STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell, char *camera)
    109108{
    110109   int k;
    … …  
    124123/****************************************************************/
    125124/* burn_apply(): Subtract the trail fits from the image */
    126 STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell)
     125STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell, char *camera)
    127126{
    128127   int k;
    … …  
    143142/* burn_check(): Find all the burned patches in a cell */
    144143STATIC int burn_check(int nx, int ny, int NX, int NY, DTYPE *data,
    145                        MTYPE *mask, CELL *cell)
     144                       MTYPE *mask, CELL *cell, char *camera)
    146145{
    147146   int err,k;
    … …  
    153152         continue;
    154153      }
    155       err = burn_test(nx, ny, NX, data, cell->rms, mask, cell->burn+k);
     154      err = burn_test(nx, ny, NX, data, cell->rms, mask, cell->burn+k, camera);
    156155   }
    157156
    … …  
    240239 */
    241240STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms,
    242                        MTYPE *mask, OBJBOX *box)
     241                       MTYPE *mask, OBJBOX *box, char *camera)
    243242{
    244243   int i, j, nburn, nref, nmed;
    245    int y0, y1, ymid, xburn, xref, dx=0, dy=0;
     244   int y0, y1, x0, x1, xwidth, ymid, xburn, xref, dx=0, dy=0;
    246245
    247246   /*introduce a maximum to the burn contrast check, to avoid a very large number in crowded regions*/
    … …  
    275274      y0 = ymid + MAX(box->ey-ymid, FIT_EDGE);
    276275      y1 = MIN(y0+100, ny-1);
     276      x0 = box->sx;
     277      x1 = box->ex;
    277278
    278279      if(2*box->ex - box->sx + 1 < nx-1) {      /* Work right? */
    … …  
    294295      dx = 0;
    295296      dy = MAX(ymid-box->sy, box->ey-ymid);
     297     
     298      if(!strcmp(camera,"gpc2")) {
     299        /*start right at the edge of the star box might not be the best*/
     300        /*add 20% of box size on top, with a maximum of ten pixels*/
     301        dy = dy + MIN((box->ey-box->sy)/5,10);
     302
     303        /*Use the inner third of the star width for this test. Should be brightest there*/
     304        xwidth = MAX((box->ex-box->sx)/3,10);     
     305        x0 = box->sx + (box->ex-box->sx)/2 - xwidth/2;
     306        x1 = box->sx + (box->ex-box->sx)/2 + xwidth/2;
     307      } else {
     308        //regular gpc1 setup
     309        x0 = box->sx;
     310        x1 = box->ex;
     311      }
     312     
    296313      if(dy < FIT_EDGE) dy = FIT_EDGE;
    297314      y0 = ymid + dy;
    298315      y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1);
    299316      if(y1-y0 > test_box) y1 = y0 + test_box;
     317     
     318   }
     319   if(VERBOSE & VERB_BOXGROW) {
     320     printf("test input dx,dy,x0,x1, y0, y1, test_box, ymid: %8d %8d %8d %8d %8d %8d %8d %8d\n", dx,dy,x0,x1, y0, y1, test_box, ymid);
    300321   }
    301322
    … …  
    305326      xburn = xref = 0;
    306327      nburn = nref = 0;
    307       for(i=box->sx; i<=box->ex; i++) {
     328      for(i=x0; i<=x1; i++) {
    308329             if(mask[i+j*NX] == MASK_NONE) {
    309330               xburn += data[i+j*NX];
    … …  
    326347      if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref;
    327348      if(VERBOSE & VERB_BOXGROW) {
    328         printf("contrast: %8d %8d %8d %8d %8d\n", j, xburn, (2*ymid-j), xref, xburn - xref);
     349        printf("contrast: %8d %8d %8d %8d %8d %8d %8d\n", j, xburn,nburn, (2*ymid-j), xref,nref, xburn - xref);
    329350      }
    330351
    … …  
    332353   box->diff = int_median(nmed, median_buf);
    333354   box->burned = (box->diff > burn_contrast);
     355   if(VERBOSE & VERB_BOXGROW) {
     356     printf("test output: %8d %8d %8d\n", box->burned, burn_contrast, box->diff);
     357   }
     358
    334359   return(0);
    335360}
    … …  
    344369/* burn_test() tests whether a box is burned or not */
    345370STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms,
    346                        MTYPE *mask, OBJBOX *box)
     371                       MTYPE *mask, OBJBOX *box, char *camera)
    347372{
    348373   int i, j, n, nrow, nmed;

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c

    r42392 r42753  
    2929   HeaderUnit ehu;
    3030   const char* ifilename = "-";
    31    char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1];
     31   char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1], *camera="gpc1";
    3232   int cellmask[MAXCELL];
    3333   int otanum;
    … …  
    176176      } else if(strncmp(argv[i], "infits=", 7) == 0) { /* infits=fname */
    177177         persistfitsfile = argv[i] + 7;
     178         
     179/* Input camera keyword (optional) */
     180      } else if(strncmp(argv[i], "camera=", 7) == 0) {  /* camera=fname */
     181         camera = argv[i] + 7;
    178182
    179183/* Keep persistence streaks which had a bad slope? */
    … …  
    577581/* Use the table-driven fits instead of calculating new ones */
    578582            burn_apply(naxis1-ovrscan1, naxis2-ovrscan2, naxis1,
    579                        buf, OTA+cell);
     583                       buf, OTA+cell,camera);
    580584/* Tell us about it? */
    581585            if(VERBOSE & VERB_NORM) burn_blab(OTA+cell);
    … …  
    584588/* Restore the old burns */
    585589            burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1,
    586                          buf, OTA+cell);
     590                         buf, OTA+cell,camera);
    587591
    588592         } else {
    589593/* Fix up the burns */
    590594            burn_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, naxis2, buf,
    591                      OTA+cell, cell);
     595                     OTA+cell, cell,camera);
    592596
    593597/* Collect up a good star list */
    … …  
    602606
    603607            persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf,
    604                            OTA+cell);
     608                           OTA+cell,camera);
    605609/* Tell us about it? */
    606610            if(VERBOSE & VERB_NORM) persist_blab(OTA+cell);

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h

    r42392 r42753  
    172172STATIC int mem_init(int nx, int ny, int NX, int NY);
    173173STATIC int burn_fix(int nx, int ny, int stride, int NY, IMTYPE *buf,
    174                     CELL *cell, int cellnum);
     174                    CELL *cell, int cellnum,char *camera);
    175175STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms,
    176                      MTYPE *mask, OBJBOX *box);
    177 STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell);
    178 STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell);
     176                     MTYPE *mask, OBJBOX *box,char *camera);
     177STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell,char *camera);
     178STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell,char *camera);
    179179STATIC int persist_read(CELL *cell, const char *infile, int apply);
    180180STATIC int persist_write(CELL *cell, const char *outfile);
    181 STATIC int persist_fix(int nx, int ny, int stride, IMTYPE *buf, CELL *cell);
     181STATIC int persist_fix(int nx, int ny, int stride, IMTYPE *buf, CELL *cell,char *camera);
    182182STATIC int persist_merge(CELL *cell);
    183183
    … …  
    187187
    188188STATIC int star_detect(int nx, int ny, int NX, int NY, DTYPE *data,
    189                        MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum);
     189                       MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum,char *camera);
    190190STATIC int burn_check(int nx, int ny, int stride, int NY, DTYPE *buf,
    191                        MTYPE *mask, CELL *cell);
     191                       MTYPE *mask, CELL *cell,char *camera);
    192192
    193193STATIC int cell_stats(int nx, int ny, int NX, int NY, DTYPE *data, CELL *cell);
    … …  
    197197
    198198STATIC int fit_trail(int nx, int ny, int NX, DTYPE *data, MTYPE *mask,
    199                     OBJBOX *box, int up, int sky, int rms, int fitfunc);
     199                    OBJBOX *box, int up, int sky, int rms, int fitfunc,char *camera);
    200200STATIC int sub_fit(int nx, int ny, int NX, IMTYPE *buf, OBJBOX *box, int sign);
    201201

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c

    r42392 r42753  
    1818/* burn_fix must have been run first! */
    1919STATIC int persist_fix(int nx, int ny, int NX, IMTYPE *buf,
    20                     CELL *cell)
     20                    CELL *cell, char *camera)
    2121{
    2222   int k, err=0;
    … …  
    4343/* Fit the trail */
    4444     err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0,
    45              cell->sky+cell->bias, cell->rms, BURN_EXP);
     45             cell->sky+cell->bias, cell->rms, BURN_EXP, camera);
    4646
    4747/* Subtract out the fit */

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c

    r42392 r42753  
    1717/* star_detect(): Find all the stars and burned patches in a cell */
    1818STATIC int star_detect(int nx, int ny, int NX, int NY, DTYPE *data,
    19                        MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum)
     19                       MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum, char *camera)
    2020{
    2121   int i, j, k, l;
    … …  
    124124/* But take a hard look at stars with really bright centers... */
    125125                if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) {
    126                  burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox);
     126                 burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox, camera);
    127127                }
    128128                boxbuf[nbox].time = cell->time;

  • branches/eam_branches/ipp-20230313/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c

    r42392 r42753  
    1414#include "burnparams.h"
    1515
    16 STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE];
     16STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE],xbuf[MAXSIZE];
    1717
    1818/****************************************************************/
    … …  
    9292STATIC int fit_trail(int nx, int ny, int NX,
    9393                    DTYPE *data, MTYPE *mask, OBJBOX *box,
    94                     int up, int bckgnd, int rms, int fitfunc)
     94                    int up, int bckgnd, int rms, int fitfunc,char *camera)
    9595{
    96    int i, j, k, err, nfit,firstfitpix;
    97    int xs, xe, y0, y1, y2, dy, yfit;
    98    double slope, zero, zsum, wsum, trial,zfun;
     96   int i, j, k, err, nfit,firstfitpix,secondfitpix, xwidth;
     97   int xs, xe, y0, y1, y2, dy, yfit,x0,x1;
     98   double slope, zero, zsum, wsum, trial,zfun,xsum,xmean,xstd;
    9999   char fooname[60];
    100100   FILE *fp=stdout;
    … …  
    104104      return(-1);
    105105   }
    106 
     106   
    107107   /*xs = box->sxfit; */
    108108   /*xe = box->exfit; */
    … …  
    110110   xs = box->sx;
    111111   xe = box->ex;
     112
    112113   
    113114   if(up) {
    … …  
    119120      yfit = box->sy;      /* Start of validity, TBD later by col */
    120121      firstfitpix = MAX((box->ey - box->sy)*2,20) ;  /* Number of pixels to include in the initial fit*/
     122      secondfitpix = 200;  /* Number of pixels to include in the second fit*/
    121123   } else {
    122124      //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
    … …  
    127129      yfit = box->ey;      /* Start of validity, TBD later by col */
    128130      firstfitpix = 100;  /* Number of pixels to include in the initial fit*/
    129 
     131      secondfitpix = 200;  /* Number of pixels to include in the second fit*/
     132   }
     133
     134   if(!strcmp(camera,"gpc2")) {
     135     /*Use the inner third of the star width for this test. Should be brightest there*/
     136     xwidth = MAX((box->ex-box->sx)/3,10);     
     137     x0 = box->sx + (box->ex-box->sx)/2 - xwidth/2;
     138     x1 = box->sx + (box->ex-box->sx)/2 + xwidth/2;
     139     /*when fitting a burn, do not want to include part of the star*/
     140     /*add 20% of box size on top, with a maximum of ten pixels*/
     141     if(up) y1 = box->ey + MIN((box->ey-box->sy)/5,10);
     142   } else {
     143     //regular gpc1 setup
     144     x0 = box->sx;
     145     x1 = box->ex;
    130146   }
    131147
    132148   if(VERBOSE & VERB_FIT) {
    133       printf("Fit requested xs,xe,y0,y1,y2,dy,up,f= %d %d %d %d %d %d %d %d\n",
    134              xs, xe, y0, y1, y2, dy, up, fitfunc);
     149      printf("Fit requested xs,xe,y0,y1,y2,dy,up,f,x0,x1= %d %d %d %d %d %d %d %d %d %d\n",
     150             xs, xe, y0, y1, y2, dy, up, fitfunc,x0,x1);
    135151   }
    136152
    … …  
    184200      ybuf[j] = Y_SCALE * dy * (j-y0);
    185201      zbuf[j]=0.0;
    186       for(i=xs, k=0; i<=xe; i++) {
     202      for(i=x0, k=0; i<=x1; i++) {
    187203             if(data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO) ) {
    188204              zbuf[j] += data[i+j*NX];
    … …  
    221237/* Trim away any really big positive deviations */
    222238   for(j=y1; j<=y2; j++) {
    223 //      if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV &&
    224 //       zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
    225239      if(abs(zbuf[j]-(ybuf[j]*slope+zero)) > MAXDEV) wbuf[j] = 0;
    226       if(VERBOSE & VERB_FITPROF) {
    227          fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j);
     240      if(j < (y2-secondfitpix) && !up) wbuf[j] = 0;
     241      if(j > (y1+secondfitpix) && up) wbuf[j] = 0;
     242      if(VERBOSE & VERB_FIT) {
     243         printf("trim j : %d %9.4f %9.4f %9.4f %d %9.4f %9.4f\n", j,ybuf[j],zbuf[j],(ybuf[j]*slope+zero),abs(zbuf[j]-(ybuf[j]*slope+zero)), wbuf[j], wbuf0[j]);
    228244      }
    229245   }
    … …  
    237253
    238254   if(VERBOSE & VERB_FITPROF) {
    239       printf("jt %d  %.3f %.3f\n", xs, zero, slope);
     255      printf("2nd fit nfit, zero, slope: %d  %.3f %.3f\n", nfit, zero, slope);
    240256      fclose(fp);
    241257   }
    … …  
    249265      linearrms(nfit, ybuf+y1, zbuf+y1, slope, zero, &trial);
    250266/* 100203 JT: bad idea: appears to be a bug in read/writing ABS(nfit) */
     267
    251268//      if(trial > 2*rms) box->nfit = -box->nfit;
    252269      if(trial > 2*rms) box->func = BURN_POSSLOPE;
    … …  
    256273   }
    257274
    258 /* Reconstruct the fit for comparison with the data col by col */
    259    for(j=y1; j<=y2; j++) zbuf[j] = exp(ybuf[j]*slope);
    260    /* And tag on the extras to determine the end of the fit */
    261    for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
    262       if(fitfunc == BURN_EXP) {
    263              zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
     275/* find where the fitted function becomes less than half the rms and cut the trail there*/
     276   for(j=(up?y1:y2); dy*j<dy*(up?y2:y1); j+=dy) {
     277   //for(j=y1; j<=y2; j++) {
     278      zfun = exp((ybuf[j]*slope)+zero);
     279      //printf("y zfun thresh : %d %9.4f %9.4f %9.4f\n",j,zfun,NEGLIGIBLE_TRAIL*rms,NEGLIGIBLE_TRAIL*rms*0.5);
     280      if(up) {
     281        if(zfun < NEGLIGIBLE_TRAIL*rms*0.5) break;
    264282      } else {
    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;
     283        if(zfun < NEGLIGIBLE_TRAIL*rms) break;     
     284      }
    274285   }   
    275286   //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 */
    279    box->slope = slope;
    280    box->nfit = 0;
    281 /* Save individual scaled versions for each column */
    282    for(i=xs; i<=xe; i++) {
    283       zsum = wsum = 0.0;
    284       for(j=y1; j<=y2; j++) {
     287   //if(j < y2 && !up) box->eyfit = j;
     288   box->eyfit = j;
     289
     290
     291   if(!strcmp(camera,"gpc2")) {
     292
     293     /*Conpute the mean and standard deviation of the trail along x*/
     294     for(i=xs; i<=xe; i++) {
     295      xbuf[i] = 0.0;
     296      for(j=y1, k=0; j<=y2; j++) {
     297        if(wbuf[j] > 0 && data[i+j*NX] != NODATA &&
     298        (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) {
     299          xbuf[i] += data[i+j*NX];
     300          k++;
     301        }
     302      }
     303      xbuf[i] = xbuf[i] / MAX(1,k);
     304     }
     305   
     306     xsum = wsum = 0.;
     307     for(i=xs; i<=xe; i++) {
     308      xsum += i*xbuf[i];
     309      wsum += xbuf[i];
     310     }
     311     xmean = xsum / wsum;
     312
     313     xsum = wsum = 0.;
     314     for(i=xs; i<=xe; i++) {
     315      xsum += xbuf[i]*pow(i-xmean,2);
     316      wsum += xbuf[i];
     317     }
     318     xstd = sqrt(xsum / (wsum-1.));
     319     if(VERBOSE & VERB_FIT) {
     320       printf("x mean xstd : %9.4f %9.4f %9.4f %9.4f\n",  xmean, xstd,xsum,wsum); 
     321     }
     322
     323     /* Update fit ranges */
     324     if(xstd > 1 && xstd < 50) {
     325      box->sxfit = round(xmean-(xstd*1.1));
     326      box->exfit = round(xmean+(xstd*1.1));
     327     } else {
     328      box->sxfit = xs;
     329      box->exfit = xs - 1;
     330      box->fiterr = FIT_ALL_GONE;
     331     }
     332   } else {
     333     /* Reconstruct the fit for comparison with the data col by col */
     334     for(j=y1; j<=y2; j++) zbuf[j] = exp((ybuf[j]*slope));
     335     /* And tag on the extras to determine the end of the fit */
     336     for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
     337      if(fitfunc == BURN_EXP) {
     338             zbuf[j] = exp((Y_SCALE*dy*(j-y0)*slope));
     339      } else {
     340             if(dy*(j-y0) <= 0) zbuf[j] = exp((log(Y_SCALE*1.0)*slope));
     341             else zbuf[j] = exp((log(Y_SCALE*dy*(j-y0))*slope));
     342      }
     343     }
     344
     345     /* Determine good x-ranges and Save the fit information */
     346     box->slope = slope;
     347     box->nfit = 0;
     348     /* Save individual scaled versions for each column */
     349     for(i=xs; i<=xe; i++) {
     350       zsum = wsum = 0.0;
     351       for(j=y1; j<=y2; j++) {
    285352             if(wbuf[j] > 0 && data[i+j*NX] != NODATA &&
    286353             (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];
    290354               wsum += zbuf[j];
    291355               zsum += data[i+j*NX] - bckgnd;
    292356             }
    293       }
    294       if(zsum < 0 || wsum <= 0) {
     357       }
     358       if(zsum < 0 || wsum <= 0) {
    295359             zsum = 0.0;
    296       } else {
     360       } else {
    297361             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 */
    303       if(zsum > NEGLIGIBLE_TRAIL*rms ||
     362       }
     363       /* FIXME: what's a really good criterion for negligible fit? */
     364       /* 100 pixels up fit is zsum or ~zsum/e */
     365       /* 100113: but also evaluate at the end of the fit for stubby trails */
     366       //printf("scaled fit comparison x zsum zbuf[y2]*zsum thresh: %d %9.4f %9.4f %9.4f %f %d %f\n", i, zsum, zbuf[y2], zbuf[y2]*zsum, NEGLIGIBLE_TRAIL, rms, NEGLIGIBLE_TRAIL*rms); 
     367
     368       if(zsum > NEGLIGIBLE_TRAIL*rms ||
    304369           zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
    305370        /* Ascertain the starting point of where the fit is good */
    … …  
    314379             box->zero[box->nfit] = zsum;
    315380             box->xfit[box->nfit] = i;
     381             //printf("y comparison : %d %d %9.4f %9.4f %d %d\n", i, box->yfit[box->nfit], box->zero[box->nfit],zsum, box->xfit[box->nfit],box->nfit); 
    316382             box->nfit += 1;
    317       }
    318    }
    319 /* Update fit ranges */
    320    if(box->nfit > 0) {
    321       box->sxfit = box->xfit[0];
    322       box->exfit = box->xfit[box->nfit-1];
    323    } else {
    324       box->sxfit = xs;
    325       box->exfit = xs - 1;
    326       box->fiterr = FIT_ALL_GONE;
     383       }
     384     }
     385     /* Update fit ranges */
     386     if(box->nfit > 0) {
     387        box->sxfit = box->xfit[0];
     388        box->exfit = box->xfit[box->nfit-1];
     389     } else {
     390        box->sxfit = xs;
     391        box->exfit = xs - 1;
     392        box->fiterr = FIT_ALL_GONE;
     393     }
     394   }
     395
     396
     397
     398   if(VERBOSE & VERB_FIT) {
     399      printf("Final fit params xsfit,xefit,y1,y2, yfit,slope,nfit = %d %d %d %d %d %.3f %d\n",
     400             box->sxfit, box->exfit, y1, y2,box->eyfit, box->slope, box->nfit);
    327401   }
    328402
Note: See TracChangeset for help on using the changeset viewer.