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