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) $
8 **********************************************************************
10 THIS FUNCTION IS USED IN DIVERGENT BEAM CONVOLUTION ALGORITHM
11 TO CALCULATE THE CONVOLUTING FUNCTIONS Q1(M) AND Q2(M) FOR
28 BOOLEAN dcon_class::dconft(REAL* q1, REAL* q2, INTEGER nconrp, REAL alfa, INTEGER nf)
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');
79 // INITIALIZE CONSTANTS FOR CONVOLUTING FUNCTION
81 for(m = 0; m < nconrp; m++)
87 // FILTER FUNCTION DEFINED BY LINES WHICH CONNECTING THE COORDINATES
88 // OF (AA,YA) AND (BB,YB)
90 if(nf != hline) goto L50;
91 fprintf(output, "\n filter function composed of following lines");
95 bb = InFile.getnum(TRUE, &eol);
98 yb = InFile.getnum(FALSE, &eol);
100 if(eol) yb = InFile.getnum(TRUE, &eol);
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);
114 if((ya < yb) || (bb > 1) || (yb < 0)) goto L25;
116 a = aa / (REAL) 2.0 / alfa;
117 b = bb / (REAL) 2.0 / alfa;
118 slope = (yb - ya) / (b - a);
121 api = aa * Consts.pi;
122 bpi = bb * Consts.pi;
125 fprintf(output, "\n (%6.4f,%6.4f) to (%6.4f,%6.4f)", aa, ya, bb, yb);
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)
131 for(m = 1; m < nconrp; m++)
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);
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;
152 bb = InFile.getnum(FALSE, &eol);
154 if(eol) bb = InFile.getnum(TRUE, &eol);
157 yb = InFile.getnum(FALSE, &eol);
159 if(eol) yb = InFile.getnum(TRUE, &eol);
166 // FILTER BANDWIDTH - A=C/ALFA
169 c = InFile.getnum(TRUE, &eol);
171 if((c <= 0.) || (c > 1.)) c = 1.;
177 //.... HAMMING LOW PASS FILTER
178 //.... F(X)=ALPHA+(1.-ALPHA)*COS(2.*PI*X/A)
180 if(nf != hhamm) goto L90;
181 alpha = InFile.getnum(FALSE, &eol);
184 fprintf(output, "\n hamming low pass filter");
185 fprintf(output, "\n alpha= %7.3f", alpha);
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)
194 for(m = is; m < nconrp; m++) //(JD 11/06/03)
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);
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;
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;
221 if(nf != hpara) goto L110;
222 fprintf(output, "\n parabolic low pass filter");
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)
227 for(m = 1; m < nconrp; m++) {
228 scmp = (REAL) sin(cmp);
229 ccmp = (REAL) cos(cmp);
230 b = (REAL) sin(sigma);
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;
241 if(nf != hexpo) goto L130;
242 fprintf(output, "\n exponential low pass filter");
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++)
249 scmp = (REAL) sin(cmp);
250 ccmp = (REAL) cos(cmp);
251 b = (REAL) sin(sigma);
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);
262 if(nf != hsinc) goto L160;
263 fprintf(output, "\n sinc low pass filter");
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)
271 for(m = is; m < nconrp; m++)
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;
285 q2[m] = (REAL) 0.5 * c * cin(Consts.pi + Consts.pi) / b / alfa;
294 if(nf != hcosi) goto L200;
295 fprintf(output, "\n cosine low pass filter");
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)
305 for(m=is; m < nconrp; m++) //(JD 11/06/03)
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;
322 q1[m] = (REAL) -0.5 / b / b;
323 q2[m] = (REAL) 0.125 * ((REAL) 4.0 + Consts.pisq) * c / alfa / b;
330 // SHEPP AND LOGAN'S SINC FILTER
333 if(nf != hshep) goto L240;
334 fprintf(output, "\n shepp and logan's sinc filter");
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)
339 is = 1; //(JD 11/06/03)
343 for(m = is; m < nconrp; m++) //(JD 11/06/03)
345 b = (REAL) sin(sigma);
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));
361 q2[m] = c / alfa / b * (cin(Consts.pi) + (REAL) 1.0);
370 // UNKNOWN FILTER NAME
371 fprintf(output, "\n can not recognize filter name -- %s job abolished", int2str(nf));
377 fprintf(output, "\n *****insufficient data items - job aborted*****");
383 fprintf(output, "\n band-width = %7.3f * nyquist cutoff freq.", c);