Changeset 42724 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c

Timestamp:

Oct 15, 2024, 10:57:32 AM (21 months ago)

Author:

tdeboer

Message:

updated working for gpc2 burntool, inclduing adding camera keyword

File:

• trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c (modified) (12 diffs)

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

@@ -14 +14 @@
 #include "burnparams.h"
 
-STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE];
+STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE],xbuf[MAXSIZE];
 
 /****************************************************************/
@@ -92 +92 @@
 STATIC int fit_trail(int nx, int ny, int NX, 
                     DTYPE *data, MTYPE *mask, OBJBOX *box,
-                    int up, int bckgnd, int rms, int fitfunc)
+                    int up, int bckgnd, int rms, int fitfunc,char *camera)
 {
-   int i, j, k, err, nfit,firstfitpix;
-   int xs, xe, y0, y1, y2, dy, yfit;
-   double slope, zero, zsum, wsum, trial,zfun;
+   int i, j, k, err, nfit,firstfitpix,secondfitpix, xwidth;
+   int xs, xe, y0, y1, y2, dy, yfit,x0,x1;
+   double slope, zero, zsum, wsum, trial,zfun,xsum,xmean,xstd;
    char fooname[60];
    FILE *fp=stdout;
@@ -104 +104 @@
       return(-1);
    }
-
+   
    /*xs = box->sxfit; */
    /*xe = box->exfit; */
@@ -110 +110 @@
    xs = box->sx;
    xe = box->ex;
+
    
    if(up) {
@@ -119 +120 @@
       yfit = box->sy;      /* Start of validity, TBD later by col */
       firstfitpix = MAX((box->ey - box->sy)*2,20) ;  /* Number of pixels to include in the initial fit*/
+      secondfitpix = 200;  /* Number of pixels to include in the second fit*/
    } else {
       //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
@@ -127 +129 @@
       yfit = box->ey;      /* Start of validity, TBD later by col */
       firstfitpix = 100;  /* Number of pixels to include in the initial fit*/
-
+      secondfitpix = 200;  /* Number of pixels to include in the second fit*/
+   }
+
+   if(!strcmp(camera,"gpc2")) {
+     /*Use the inner third of the star width for this test. Should be brightest there*/
+     xwidth = MAX((box->ex-box->sx)/3,10);      
+     x0 = box->sx + (box->ex-box->sx)/2 - xwidth/2;
+     x1 = box->sx + (box->ex-box->sx)/2 + xwidth/2;
+     /*when fitting a burn, do not want to include part of the star*/
+     /*add 20% of box size on top, with a maximum of ten pixels*/
+     if(up) y1 = box->ey + MIN((box->ey-box->sy)/5,10);
+   } else {
+     //regular gpc1 setup
+     x0 = box->sx;
+     x1 = box->ex;
    }
 
    if(VERBOSE & VERB_FIT) {
-      printf("Fit requested xs,xe,y0,y1,y2,dy,up,f= %d %d %d %d %d %d %d %d\n", 
-             xs, xe, y0, y1, y2, dy, up, fitfunc);
+      printf("Fit requested xs,xe,y0,y1,y2,dy,up,f,x0,x1= %d %d %d %d %d %d %d %d %d %d\n", 
+             xs, xe, y0, y1, y2, dy, up, fitfunc,x0,x1);
    }
 
@@ -184 +200 @@
       ybuf[j] = Y_SCALE * dy * (j-y0);
       zbuf[j]=0.0;
-      for(i=xs, k=0; i<=xe; i++) {
+      for(i=x0, k=0; i<=x1; i++) {
              if(data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO) ) {
               zbuf[j] += data[i+j*NX];
@@ -221 +237 @@
 /* Trim away any really big positive deviations */
    for(j=y1; j<=y2; j++) {
-//      if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV && 
-//       zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
       if(abs(zbuf[j]-(ybuf[j]*slope+zero)) > MAXDEV) wbuf[j] = 0;
-      if(VERBOSE & VERB_FITPROF) {
-         fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j);
+      if(j < (y2-secondfitpix) && !up) wbuf[j] = 0;
+      if(j > (y1+secondfitpix) && up) wbuf[j] = 0;
+      if(VERBOSE & VERB_FIT) {
+         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]);
       }
    }
@@ -237 +253 @@
 
    if(VERBOSE & VERB_FITPROF) {
-      printf("jt %d  %.3f %.3f\n", xs, zero, slope);
+      printf("2nd fit nfit, zero, slope: %d  %.3f %.3f\n", nfit, zero, slope);
       fclose(fp);
    }
@@ -249 +265 @@
       linearrms(nfit, ybuf+y1, zbuf+y1, slope, zero, &trial);
 /* 100203 JT: bad idea: appears to be a bug in read/writing ABS(nfit) */
+
 //      if(trial > 2*rms) box->nfit = -box->nfit;
       if(trial > 2*rms) box->func = BURN_POSSLOPE;
@@ -256 +273 @@
    }
 
-/* Reconstruct the fit for comparison with the data col by col */
-   for(j=y1; j<=y2; j++) zbuf[j] = exp(ybuf[j]*slope);
-   /* And tag on the extras to determine the end of the fit */
-   for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
-      if(fitfunc == BURN_EXP) {
-             zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
+/* find where the fitted function becomes less than half the rms and cut the trail there*/
+   for(j=(up?y1:y2); dy*j<dy*(up?y2:y1); j+=dy) {
+   //for(j=y1; j<=y2; j++) {
+      zfun = exp((ybuf[j]*slope)+zero);
+      //printf("y zfun thresh : %d %9.4f %9.4f %9.4f\n",j,zfun,NEGLIGIBLE_TRAIL*rms,NEGLIGIBLE_TRAIL*rms*0.5);
+      if(up) {
+        if(zfun < NEGLIGIBLE_TRAIL*rms*0.5) break;
       } else {
-             if(dy*(j-y0) <= 0) zbuf[j] = exp(log(Y_SCALE*1.0)*slope);
-             else zbuf[j] = exp(log(Y_SCALE*dy*(j-y0))*slope);
-      }
-   }
-
-/* find where the fitted function becomes less than half the rms and cut the trail there*/
-   for(j=y1; j<=y2; j++) {
-      zfun = exp((ybuf[j]*slope)+zero);
-      if(zfun > NEGLIGIBLE_TRAIL*rms) break;
+        if(zfun < NEGLIGIBLE_TRAIL*rms) break;      
+      }
    }   
    //For burns, just let it go to the top
-   if(j < y2 && !up) box->eyfit = j;
-
-/* Determine good x-ranges and Save the fit information */
-   box->slope = slope;
-   box->nfit = 0;
-/* Save individual scaled versions for each column */
-   for(i=xs; i<=xe; i++) {
-      zsum = wsum = 0.0;
-      for(j=y1; j<=y2; j++) {
+   //if(j < y2 && !up) box->eyfit = j;
+   box->eyfit = j;
+
+
+   if(!strcmp(camera,"gpc2")) {
+
+     /*Conpute the mean and standard deviation of the trail along x*/
+     for(i=xs; i<=xe; i++) {
+      xbuf[i] = 0.0;
+      for(j=y1, k=0; j<=y2; j++) {
+        if(wbuf[j] > 0 && data[i+j*NX] != NODATA && 
+        (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) {
+          xbuf[i] += data[i+j*NX];
+          k++;
+        }
+      }
+      xbuf[i] = xbuf[i] / MAX(1,k);
+     }
+    
+     xsum = wsum = 0.;
+     for(i=xs; i<=xe; i++) {
+      xsum += i*xbuf[i];
+      wsum += xbuf[i];
+     }
+     xmean = xsum / wsum;
+
+     xsum = wsum = 0.;
+     for(i=xs; i<=xe; i++) {
+      xsum += xbuf[i]*pow(i-xmean,2);
+      wsum += xbuf[i];
+     }
+     xstd = sqrt(xsum / (wsum-1.));
+     if(VERBOSE & VERB_FIT) {
+       printf("x mean xstd : %9.4f %9.4f %9.4f %9.4f\n",  xmean, xstd,xsum,wsum);  
+     }
+
+     /* Update fit ranges */
+     if(xstd > 1 && xstd < 50) {
+      box->sxfit = round(xmean-(xstd*1.1));
+      box->exfit = round(xmean+(xstd*1.1));
+     } else {
+      box->sxfit = xs;
+      box->exfit = xs - 1;
+      box->fiterr = FIT_ALL_GONE;
+     }
+   } else {
+     /* Reconstruct the fit for comparison with the data col by col */
+     for(j=y1; j<=y2; j++) zbuf[j] = exp((ybuf[j]*slope));
+     /* And tag on the extras to determine the end of the fit */
+     for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
+      if(fitfunc == BURN_EXP) {
+             zbuf[j] = exp((Y_SCALE*dy*(j-y0)*slope));
+      } else {
+             if(dy*(j-y0) <= 0) zbuf[j] = exp((log(Y_SCALE*1.0)*slope));
+             else zbuf[j] = exp((log(Y_SCALE*dy*(j-y0))*slope));
+      }
+     }
+
+     /* Determine good x-ranges and Save the fit information */
+     box->slope = slope;
+     box->nfit = 0;
+     /* Save individual scaled versions for each column */
+     for(i=xs; i<=xe; i++) {
+       zsum = wsum = 0.0;
+       for(j=y1; j<=y2; j++) {
              if(wbuf[j] > 0 && data[i+j*NX] != NODATA && 
              (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) {
-//             if(wbuf[j] > 0 && data[i+j*NX] != NODATA) {
-//             wsum += wbuf[j];
-//             zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j];
                wsum += zbuf[j];
                zsum += data[i+j*NX] - bckgnd;
              }
-      }
-      if(zsum < 0 || wsum <= 0) {
+       }
+       if(zsum < 0 || wsum <= 0) {
              zsum = 0.0;
-      } else {
+       } else {
              zsum = zsum / wsum;
-      }
-      /* FIXME: what's a really good criterion for negligible fit? */
-      /* 100 pixels up fit is zsum or ~zsum/e */
-//    if(zsum > NEGLIGIBLE_TRAIL*rms) {
-      /* 100113: but also evaluate at the end of the fit for stubby trails */
-      if(zsum > NEGLIGIBLE_TRAIL*rms ||
+       }
+       /* FIXME: what's a really good criterion for negligible fit? */
+       /* 100 pixels up fit is zsum or ~zsum/e */
+       /* 100113: but also evaluate at the end of the fit for stubby trails */
+       //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);  
+
+       if(zsum > NEGLIGIBLE_TRAIL*rms ||
            zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
         /* Ascertain the starting point of where the fit is good */
@@ -314 +379 @@
              box->zero[box->nfit] = zsum;
              box->xfit[box->nfit] = i;
+             //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);  
              box->nfit += 1;
-      }
-   }
-/* Update fit ranges */
-   if(box->nfit > 0) {
-      box->sxfit = box->xfit[0];
-      box->exfit = box->xfit[box->nfit-1];
-   } else {
-      box->sxfit = xs;
-      box->exfit = xs - 1;
-      box->fiterr = FIT_ALL_GONE;
+       }
+     }
+     /* Update fit ranges */
+     if(box->nfit > 0) {
+        box->sxfit = box->xfit[0];
+        box->exfit = box->xfit[box->nfit-1];
+     } else {
+        box->sxfit = xs;
+        box->exfit = xs - 1;
+        box->fiterr = FIT_ALL_GONE;
+     }
+   }
+
+
+
+   if(VERBOSE & VERB_FIT) {
+      printf("Final fit params xsfit,xefit,y1,y2, yfit,slope,nfit = %d %d %d %d %d %.3f %d\n", 
+             box->sxfit, box->exfit, y1, y2,box->eyfit, box->slope, box->nfit);
    }