Changeset 42392 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c
- Timestamp:
- Feb 14, 2023, 12:48:54 PM (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c
r26766 r42392 14 14 #include "burnparams.h" 15 15 16 STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE] ;16 STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE]; 17 17 18 18 /****************************************************************/ … … 86 86 #define MAXWGT 2.0 /* Maximum allowed weight */ 87 87 #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 */ 89 89 #define MAXRMS 4.0 /* Minimum factor of RMS to be chopped */ 90 90 … … 94 94 int up, int bckgnd, int rms, int fitfunc) 95 95 { 96 int i, j, k, err, nfit ;96 int i, j, k, err, nfit,firstfitpix; 97 97 int xs, xe, y0, y1, y2, dy, yfit; 98 double slope, zero, zsum, wsum, trial ;98 double slope, zero, zsum, wsum, trial,zfun; 99 99 char fooname[60]; 100 100 FILE *fp=stdout; … … 105 105 } 106 106 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 110 113 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 */ 116 118 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*/ 118 121 } 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 */ 124 126 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 126 130 } 127 131 … … 147 151 box->up = up; 148 152 box->y0 = y0; 153 box->eyfit = up?y2:y1; 149 154 for(i=0; i<xe-xs+1; i++) { 150 155 box->zero[i] = 0.0; … … 157 162 box->slope = -1.0; 158 163 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 */ 160 165 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 */ 163 168 } 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; 166 171 } 167 172 return(0); … … 174 179 175 180 /* Accumulate points to fit */ 181 /* use wbuf0 to select only the first FIRSTFITPIX pixels for first fit */ 176 182 nfit = y2 - y1 + 1; 177 183 for(j=y1; j<=y2; j++) { … … 179 185 zbuf[j]=0.0; 180 186 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 } 189 194 zbuf[j] = zbuf[j] / MAX(1,k) - bckgnd; 190 /* Linearize for fit... */ 195 196 /* Linearize for fit... */ 191 197 if(fitfunc == BURN_PWR) ybuf[j] = log(ybuf[j]); 192 198 wbuf[j] = 0; 199 wbuf0[j] = 0; 193 200 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]); 197 208 } 198 209 } 199 210 /* 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); 201 212 if(VERBOSE & VERB_FIT) { 202 213 printf("nfit= %d err= %d slope= %.3f zero= %.3f\n", … … 210 221 /* Trim away any really big positive deviations */ 211 222 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; 214 226 if(VERBOSE & VERB_FITPROF) { 215 227 fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j); … … 239 251 // if(trial > 2*rms) box->nfit = -box->nfit; 240 252 if(trial > 2*rms) box->func = BURN_POSSLOPE; 253 box->sxfit = xs; 254 box->exfit = xs - 1; 241 255 return(-1); 242 256 } … … 244 258 /* Reconstruct the fit for comparison with the data col by col */ 245 259 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 */ 247 261 for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) { 248 262 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); 251 264 } 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 */ 258 279 box->slope = slope; 259 280 box->nfit = 0; … … 262 283 zsum = wsum = 0.0; 263 284 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 } 272 293 } 273 294 if(zsum < 0 || wsum <= 0) { 274 zsum = 0.0;295 zsum = 0.0; 275 296 } 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 */ 281 303 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 */ 286 306 // box->yfit[box->nfit] = yfit; /* Vanilla, misses center */ 287 307 // 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; 297 317 } 298 318 } … … 306 326 box->fiterr = FIT_ALL_GONE; 307 327 } 328 308 329 // if(box->nfit <= 0) fprintf(stderr, "WHOA, got nfit = 0\n"); 309 330 return(0);
Note:
See TracChangeset
for help on using the changeset viewer.
