src/snark/dcon_dconft.cpp

   1 /*
   2 **********************************************************************
   3 $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/dcon_dconft.cpp $
   5 $LastChangedRevision: 85 $
   6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7 $Author: agulati $
   8 **********************************************************************
   9
  10  THIS FUNCTION IS USED IN DIVERGENT BEAM CONVOLUTION ALGORITHM
  11  TO CALCULATE THE CONVOLUTING FUNCTIONS Q1(M) AND Q2(M) FOR
  12  ARC DETECTOR CASE
  13 */
  14
  15 #include <cmath>
  16 #include <cstdio>
  17
  18 #include "blkdta.h"
  19 #include "region.h"
  20 #include "consts.h"
  21 #include "uiod.h"
  22 #include "int2str.h"
  23 #include "infile.h"
  24 #include "cin.h"
  25
  26 #include "dcon.h"
  27
  28 BOOLEAN dcon_class::dconft(REAL* q1, REAL* q2, INTEGER nconrp, REAL alfa, INTEGER nf)
  29 {
  30
  31   INTEGER m;
  32   REAL aa;
  33   REAL ya;
  34   REAL bb;
  35   REAL yb;
  36   REAL a;
  37   REAL b;
  38   REAL slope;
  39   REAL sigma;
  40   REAL api;
  41   REAL bpi;
  42   REAL ampi;
  43   REAL bmpi;
  44   REAL sinsgm;
  45   REAL sinamp;
  46   REAL cosamp;
  47   REAL sinbmp;
  48   REAL cosbmp;
  49   REAL c;
  50   REAL cpi;
  51   REAL cm;
  52   REAL cmp;
  53   REAL alpha;
  54   INTEGER is;
  55   REAL scmp;
  56   REAL ccmp;
  57   REAL cm2;
  58   REAL cmp2;
  59   REAL e;
  60   REAL xcin;
  61   REAL cinf;
  62
  63   REAL y0;
  64   REAL p;
  65   REAL cr;
  66
  67
  68   BOOLEAN eol;
  69
  70   static const INTEGER hline = CHAR2INT('l','i','n','e');
  71   static const INTEGER hhamm = CHAR2INT('h','a','m','m');
  72   static const INTEGER hpara = CHAR2INT('p','a','r','a');
  73   static const INTEGER hexpo = CHAR2INT('e','x','p','o');
  74   static const INTEGER hsinc = CHAR2INT('s','i','n','c');
  75   static const INTEGER hcosi = CHAR2INT('c','o','s','i');
  76   static const INTEGER hshep = CHAR2INT('s','h','e','p');
  77
  78
  79   // INITIALIZE CONSTANTS FOR CONVOLUTING FUNCTION
  80
  81   for(m = 0; m < nconrp; m++)
  82   {
  83     q1[m] = 0.;
  84     q2[m] = 0.;
  85   }
  86
  87   // FILTER FUNCTION DEFINED BY LINES WHICH CONNECTING THE COORDINATES
  88   // OF (AA,YA) AND (BB,YB)
  89
  90   if(nf != hline) goto L50;
  91   fprintf(output, "\n          filter function composed of following lines");
  92
  93   aa = 0.;
  94   ya = 1.;
  95   bb = InFile.getnum(TRUE, &eol);
  96
  97   if(eol) goto L400;
  98   yb = InFile.getnum(FALSE, &eol);
  99
 100   if(eol) yb = InFile.getnum(TRUE, &eol);
 101   if(eol) goto L400;
 102
 103 L20:
 104   if((aa-bb) > 0)
 105   {
 106 L25:
 107     fprintf(output, "\n          unrealistic parameters aa = %10.4f   bb = %10.4f", aa, bb);
 108     fprintf(output, "\n                                 ya = %10.4f   yb = %10.4f", ya, yb);
 109         return TRUE;
 110   }
 111
 112   if((aa-bb) < 0)
 113   {
 114     if((ya < yb) || (bb > 1) || (yb < 0)) goto L25;
 115
 116     a = aa / (REAL) 2.0 / alfa;
 117     b = bb / (REAL) 2.0 / alfa;
 118     slope = (yb - ya) / (b - a);
 119     y0 = ya - a * slope;
 120     sigma = alfa;
 121     api = aa * Consts.pi;
 122     bpi = bb * Consts.pi;
 123     ampi = api;
 124     bmpi = bpi;
 125     fprintf(output, "\n            (%6.4f,%6.4f) to (%6.4f,%6.4f)", aa, ya, bb, yb);
 126
 127     p = (REAL) -4.0 * Consts.pisq * (y0*(b+a)*(b-a) / (REAL) 2.0 + slope * (b*b*b-a*a*a) / (REAL) 3.0);
 128     q1[0] = q1[0] + p;                           //(JD 11/06/03)
 129     q2[0] = q2[0] - (REAL) 2.0 * p;              //(JD 11/06/03)
 130
 131     for(m = 1; m < nconrp; m++)
 132         {
 133                 sinsgm = (REAL) sin(sigma);
 134                 sinamp = (REAL) sin(ampi);
 135                 cosamp = (REAL) cos(ampi);
 136                 sinbmp = (REAL) sin(bmpi);
 137                 cosbmp = (REAL) cos(bmpi);
 138
 139                 q1[m] -= (ya * cosamp - yb * cosbmp + slope * (sinbmp - sinamp) / Consts.twopi / sigma) / sinsgm / sinsgm;
 140                 q2[m] += (Consts.twopi * (yb * b * sinbmp - ya * a * sinamp) + slope * (b * cosbmp - a * cosamp - (sinbmp - sinamp) / Consts.twopi / sigma) / sigma) / sinsgm;
 141
 142                 sigma += alfa;
 143                 ampi += api;
 144                 bmpi += bpi;
 145     }
 146   }
 147
 148   if((yb != 0.))
 149   {
 150     aa = bb;
 151     ya = yb;
 152     bb = InFile.getnum(FALSE, &eol);
 153
 154     if(eol) bb = InFile.getnum(TRUE, &eol);
 155     if(eol) goto L400;
 156
 157     yb = InFile.getnum(FALSE, &eol);
 158
 159     if(eol) yb = InFile.getnum(TRUE, &eol);
 160
 161     if(eol) goto L400;
 162     goto L20;
 163   }
 164   return FALSE;
 165
 166   // FILTER BANDWIDTH -  A=C/ALFA
 167
 168 L50:
 169   c = InFile.getnum(TRUE, &eol);
 170   if(eol) goto L400;
 171   if((c <= 0.) || (c > 1.)) c = 1.;
 172   cpi = c * Consts.pi;
 173   sigma = alfa;
 174   cm = c;
 175   cmp = cpi;
 176
 177 //.... HAMMING LOW PASS FILTER
 178 //....   F(X)=ALPHA+(1.-ALPHA)*COS(2.*PI*X/A)
 179
 180   if(nf != hhamm) goto L90;
 181   alpha = InFile.getnum(FALSE, &eol);
 182
 183   if(eol) alpha = 1.;
 184   fprintf(output, "\n          hamming low pass filter");
 185   fprintf(output, "\n           alpha= %7.3f", alpha);
 186
 187   cr = (REAL) 1.0 / c;
 188   Region.d1 = (REAL) 1.0 - alpha;
 189   q1[0] = (REAL) 0.5 * SQR(c / alfa) * ((REAL) 4.0 * Region.d1 - alpha * Consts.pisq);  //(JD 11/06/03)
 190   q2[0] = (REAL) -2.0 * q1[0];                //(JD 11/06/03)
 191   is = 1;                                     //(JD 11/06/03)
 192
 193 L60:
 194         for(m = is; m < nconrp; m++)                    //(JD 11/06/03)
 195         {
 196                 if(fabs(((float)(m)) - cr) < Consts.zero) goto L80;          //(JD 11/06/03)
 197                 scmp = (REAL) sin(cmp);
 198                 ccmp = (REAL) cos(cmp);
 199                 b = (REAL) sin(sigma);
 200                 cm2 = cm * cm;
 201                 q1[m] = -(cm2 - alpha + ((Region.d1 - alpha) * cm2 + alpha) * ccmp) / b / b / (cm2 - (REAL) 1.0);
 202                 q2[m] = c * ((cm2 - (REAL) 1.0) * (cm2 * (alpha - Region.d1) - alpha) * Consts.pi * scmp - (REAL) 2.0 * cm * Region.d1 * ((REAL) 1.0 + ccmp)) / b / (REAL) SQR(cm2 - (REAL) 1.0) / alfa;
 203                 sigma += alfa;
 204                 cmp += cpi;
 205                 cm += c;
 206         }
 207   goto L600;
 208
 209 L80:
 210   b = (REAL) sin(sigma);
 211   q1[m] = (REAL) -2.0 * alpha / b / b;
 212   q2[m] = (REAL) 0.25 * Region.d1 * c * Consts.pisq / b / alfa;
 213   is = m+1;
 214   sigma += alfa;
 215   cm += c;
 216   cmp += cpi;
 217   goto L60;
 218
 219
 220 L90:
 221   if(nf != hpara) goto L110;
 222   fprintf(output, "\n           parabolic low pass filter");
 223
 224   q1[0] = (REAL) -0.25 * SQR(cpi/alfa);           //(JD 11/06/03)
 225   q2[0] = (REAL) -2.0 * q1[0];                    //(JD 11/06/03)
 226
 227   for(m = 1; m < nconrp; m++) {
 228     scmp = (REAL) sin(cmp);
 229     ccmp = (REAL) cos(cmp);
 230     b = (REAL) sin(sigma);
 231     cmp2 = cmp * cmp;
 232     q1[m] = ((REAL) 2.0 *(cmp * scmp - (REAL) 1.0 + ccmp) - cmp2) / (REAL) SQR(b*cmp);
 233     q2[m] = (REAL) 2.0 *((REAL) 2.0*cmp*scmp+((REAL) 2.0-cmp2)*ccmp-(REAL) 2.0) / b / sigma / cmp2;
 234     sigma += alfa;
 235     cmp += cpi;
 236   }
 237   goto L600;
 238
 239
 240 L110:
 241         if(nf != hexpo) goto L130;
 242         fprintf(output, "\n           exponential low pass filter");
 243
 244         e = (REAL) 2.7182818;
 245         q1[0] = (REAL) -0.5 * (e - (REAL) 2.0) / (e - (REAL) 1.0) * SQR(cpi/alfa);       //(JD 11/06/03)
 246         q2[0] = (REAL) -2. * q1[0];                                                      //(JD 11/06/03)
 247         for(m = 1; m < nconrp; m++)
 248         {
 249                 scmp = (REAL) sin(cmp);
 250                 ccmp = (REAL) cos(cmp);
 251                 b = (REAL) sin(sigma);
 252                 cmp2 = cmp * cmp;
 253                 q1[m] = -e / (e - (REAL) 1.0) * ((REAL) 1.0 - ccmp + cmp * ((e - (REAL) 1.0) * cmp / e - scmp)) / b / b / ((REAL) 1.0 + cmp2);
 254                 q2[m] = e / (e - (REAL) 1.0) * cmp2 * ((REAL) -2.0 / e + ((REAL) 1.0 - cmp2) * ccmp + (REAL) 2.0 * cmp * scmp) / b / sigma / SQR((REAL) 1.0 + cmp2);
 255                 sigma += alfa;
 256                 cmp += cpi;
 257         }
 258   goto L600;
 259
 260
 261 L130:
 262   if(nf != hsinc) goto L160;
 263   fprintf(output, "\n           sinc low pass filter");
 264
 265   cr = (REAL) 1.0 / c;
 266   q1[0] = (REAL) -2.0 * SQR(c/alfa);                  //(JD 11/06/03)
 267   q2[0] = (REAL) -2.0 * q1[0];                        //(JD 11/06/03)
 268   is = 1;                                             //(JD 11/06/03)
 269
 270 L140:
 271         for(m = is; m < nconrp; m++)
 272         {                  //(JD 11/06/03)
 273                 b = (REAL) sin(sigma);
 274                 xcin = cin(cmp - Consts.pi) - cin(cmp + Consts.pi);
 275                 q1[m] = (REAL) 0.5 * cm * xcin / b / b;
 276                 if(fabs(((float)(m)) - cr) < Consts.zero) goto L155;       //(JD 11/19/03)
 277                 q2[m] = c * ((REAL) -0.5 * xcin - cm * ((REAL) 1.0 + (REAL) cos(cmp)) / (cm * cm - (REAL) 1.0)) / b / alfa;
 278                 cm += c;
 279                 sigma += alfa;
 280                 cmp += cpi;
 281         }
 282   goto L600;
 283
 284 L155:
 285   q2[m] = (REAL) 0.5 * c * cin(Consts.pi + Consts.pi) / b / alfa;
 286   cm += c;
 287   sigma += alfa;
 288   cmp += cpi;
 289   is = m+1;
 290   goto L140;
 291
 292
 293 L160:
 294   if(nf != hcosi) goto L200;
 295   fprintf(output, "\n           cosine low pass filter");
 296
 297   cr = (REAL) 0.5 / c;
 298   q1[0] = (REAL) 2.0 * ((REAL) 2.0 - Consts.pi) * SQR(c / alfa);     //(JD 11/06/03)
 299   q2[0] = (REAL) -2.0 * q1[0];                                       //(JD 11/06/03)
 300   is=1;                                                              //(JD 11/06/03)
 301
 302 L170:
 303
 304
 305         for(m=is; m < nconrp; m++)                                                                              //(JD 11/06/03)
 306         {
 307                 b = (REAL) sin(sigma);
 308                 if(fabs(float(m) - cr) < Consts.zero) goto L190;     //(JD 11/19/03)
 309                 scmp = (REAL) sin(cmp);
 310                 ccmp = (REAL) cos(cmp);
 311                 cm2 = (REAL) 4.0 * cm * cm;
 312                 q1[m] = ((REAL) 2.0 * cm * scmp - cm2) / b / b / (cm2 - (REAL) 1.0);
 313                 q2[m] = (REAL) 2.0 * c / alfa * ((REAL) -4.0 * cm + (cm2 + (REAL) 1.0) * scmp - cmp * (cm2 - (REAL) 1.0) * ccmp) / SQR(cm2 - (REAL) 1.0) / b;
 314                 cm += c;
 315                 sigma += alfa;
 316                 cmp += cpi;
 317         }
 318
 319   goto L600;
 320
 321 L190:
 322   q1[m] = (REAL) -0.5 / b / b;
 323   q2[m] = (REAL) 0.125 * ((REAL) 4.0 + Consts.pisq) * c / alfa / b;
 324   cm += c;
 325   sigma += alfa;
 326   cmp += cpi;
 327   is = m + 1;
 328   goto L170;
 329
 330   // SHEPP AND LOGAN'S SINC FILTER
 331
 332 L200:
 333   if(nf != hshep) goto L240;
 334   fprintf(output, "\n           shepp and logan's sinc filter");
 335
 336   q1[0] = (REAL) -4.0 * SQR(c/alfa);                //(JD 11/06/03)
 337   q2[0] = (REAL) -2.0 * q1[0];                      //(JD 11/06/03)
 338   cr = (REAL) 0.5 / c;
 339   is = 1;                                           //(JD 11/06/03)
 340
 341 L210:
 342
 343         for(m = is; m < nconrp; m++)                                     //(JD 11/06/03)
 344         {
 345                 b = (REAL) sin(sigma);
 346                 cm2 = cm * cm;
 347                 cinf = cin(cmp + Consts.pid2) - cin(cmp - Consts.pid2);
 348                 q1[m] = -cm * cinf / b / b;
 349                 //if(fabs(((float)(m - 1)) - cr) < Consts.zero) goto L230;
 350                 if(fabs(((float)(m)) - cr) < Consts.zero) goto L230;           //(JD 11/19/03)
 351                 q2[m] = c / alfa / b * (cinf + ((REAL) 2.0 * cm2 * (REAL) sin(cmp) - cm) / (cm2 - (REAL) 0.25));
 352                 sigma += alfa;
 353                 cm += c;
 354                 cmp += cpi;
 355
 356         }
 357
 358   goto L600;
 359
 360 L230:
 361   q2[m] = c / alfa / b * (cin(Consts.pi) + (REAL) 1.0);
 362   sigma += alfa;
 363   cm += c;
 364   cmp += cpi;
 365   is = m+1;
 366   goto L210;
 367
 368 L240:
 369
 370   // UNKNOWN FILTER NAME
 371   fprintf(output, "\n           can not recognize filter name -- %s  job abolished", int2str(nf));
 372   return TRUE;
 373
 374   // WRONG INPUT DATA
 375
 376 L400:
 377   fprintf(output, "\n           *****insufficient data items - job aborted*****");
 378   return TRUE;
 379
 380   // JOB FINISHED
 381
 382 L600:
 383   fprintf(output, "\n           band-width = %7.3f * nyquist cutoff freq.", c);
 384   return FALSE;
 385 }
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395