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/bwray.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 this subroutine is the WRAY for blobs
27 void Blob_class::bwray(INTEGER np, INTEGER nr, INTEGER* list, REAL* weight,
28 INTEGER* numb, REAL* snorm)
31 INTEGER ac, bc, bc1, ac1;
45 REAL ax, ay, axid, ayid;
47 REAL theta, sinth, costh;
48 REAL prfx, prlx, prfy, prly;
51 REAL ddeltain, ddeltaout;
53 REAL ydelta; //xdelta not used
66 posit(np, nr, &ax, &ay, &mx, &my);
68 if (0 >= raylen(SOT_rect, ax, ay, mx, my, 0, 0, GeoPar.nelem * GeoPar.pixsiz / 2.0, GeoPar.nelem * GeoPar.pixsiz / 2.0, 1, 0))
70 return; // bug 111, if ray(np, nr) miss all pixels in the image, return numb=0, by wei, 6/2005
73 multf = Blob.support / Blob.delta;
75 samp_support = (REAL) samp / Blob.support;
77 laf = (INTEGER) (-0.5 * ((int) (Blob.M / 2) + (int) (Blob.H / 2)));
80 * ((int) ((Blob.M - 1) / 2) + (int) ((Blob.H - 1) / 2)));
82 lbf = -1 * (INTEGER) (Blob.H / 2);
83 lbl = (INTEGER) ((Blob.H - 1) / 2);
85 lmf = -1 * (INTEGER) (Blob.M / 2);
86 lml = (INTEGER) ((Blob.M - 1) / 2);
88 prfx = -0.25 * (Blob.M - 1.0);
91 prfy = -0.25 * Consts.sqrt3 * (Blob.H - 1.0);
94 H2 = (INTEGER) (Blob.H / 2);
95 M2 = (INTEGER) (Blob.M / 2);
96 axid = ax / Blob.delta;
97 ayid = ay / Blob.delta;
102 theta = Consts.pi * 0.5;
104 theta = (REAL) atan2(my, mx);
107 theta += Consts.pi * 2; //use atan2 instead of acos for more acuracy,the range of atan2 is [-pi,pi], by wei 1/2005
109 //added 7/19, think of theta as always being betw 0 and pi
110 if (theta > Consts.pi)
112 theta = theta - Consts.pi;
114 costh = -costh; //never used anyway
119 if (fabs(my) > Consts.zero)
120 { // non-horizontal lines
123 ord = (axid - islope * ayid);
124 coef = 0.5 * (Consts.sqrt3 * islope - 1);
126 if (theta < Consts.pi / 6 || theta > 2.0 * Consts.pi / 3.0)
129 cte = multf / my; // for the "side rays"
132 fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
134 mcov = (INTEGER) ceil(multf / fcosang);
135 ddeltain = Blob.delta * fcosang;
136 ddeltaout = Blob.delta * sinth;
138 if (theta < Consts.pi / 6)
150 yup = (x - ord + ctes) * slope;
151 ylo = (x - ord - ctes) * slope;
159 x = ord - ctes + islope * yup;
162 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
167 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
171 x = ord + ctes + islope * ylo;
172 ac = (INTEGER) floor(x - ylo * Consts.isqrt3);
179 yup = (x - ord + ctes) * slope;
180 ylo = (x - ord - ctes) * slope;
188 x = ord + ctes + islope * yup;
190 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
195 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
199 x = ord - ctes + islope * yup;
200 ac = (INTEGER) ceil(x - yup * Consts.isqrt3);
205 bc = (INTEGER) (floor(icoef * (a - ord)));
207 bf = MAX0(bc - mcov, lbf);
208 bl = MIN0(bc + mcov, lbl);
211 y = 0.5 * Consts.sqrt3 * bc;
212 dc = Blob.delta * (REAL) fabs(-my * (axid - x) + mx * (ayid - y));
225 //starting the loops--------------------------------
234 if (m >= lmf && m <= lml && b <= lbl)
237 if (n < Blob.support)
240 list[*numb] = h1 * Blob.M + m1;
241 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
243 weight[*numb] = Blob.ltable[nint];
245 *snorm += weight[*numb] * weight[*numb];
266 if (m >= lmf && m <= lml && b >= lbf)
268 if (n < Blob.support)
271 list[*numb] = h1 * Blob.M + m1;
272 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
274 weight[*numb] = Blob.ltable[nint];
276 *snorm += weight[*numb] * weight[*numb];
289 dc = dc + (REAL) sign * ddeltaout;
293 if (dc > ddeltain || dc < 0.0)
296 dc = (REAL) sign * (fabs(dc) - ddeltain);
299 bf = MAX0(bc - mcov, lbf);
300 bl = MIN0(bc + mcov, lbl);
315 fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
316 mcov = (INTEGER) ceil(multf / sinth);
317 // constants independent of the first grid point
320 ddeltain = Blob.delta * sinth;
321 ddeltaout = Blob.delta * fcosang;
323 if (theta < Consts.pi / 3)
332 //----------------------------------------------
333 // parameters of the first grid point
337 ac = (INTEGER) (floor(coef * b + ord));
339 af = MAX0(ac - mcov, laf);
340 al = MIN0(ac + mcov, lal);
343 y = 0.5 * Consts.sqrt3 * b;
344 dc = Blob.delta * (REAL) fabs(-my * (x - axid) + mx * (y - ayid));
357 // starting the loop .....................
366 if (m >= lmf && m <= lml)
369 if (n < Blob.support)
371 list[*numb] = h1 * Blob.M + m1;
372 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
374 weight[*numb] = Blob.ltable[nint];
376 *snorm += weight[*numb] * weight[*numb];
395 if (m >= lmf && m <= lml)
397 if (n < Blob.support)
400 list[*numb] = h1 * Blob.M + m1;
402 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
404 weight[*numb] = Blob.ltable[nint];
406 *snorm += weight[*numb] * weight[*numb];
417 dc = dc + (REAL) sign * ddeltaout;
421 if (dc > ddeltain || dc < 0.0)
424 dc = (REAL) sign * (fabs(dc) - ddeltain);
429 af = MAX0(ac - mcov, laf);
430 al = MIN0(ac + mcov, lal);
442 } // if non- horizontal line.
445 { //if horizontal line (my == 0)
447 fcosang = cos(Consts.pi / 6.0);
449 mcov = (INTEGER) ceil(multf / fcosang);
451 // constants independent of the first grid point
452 ddeltain = Blob.delta * fcosang;
454 //parameters of the first grid point
459 ylo = yup - 2 * multf;
461 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
468 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
471 //done with af and al
474 bc = (INTEGER) (floor(2 * ayid * Consts.isqrt3));
475 bf = MAX0(bc - mcov, lbf);
476 bl = MIN0(bc + mcov, lbl);
478 y = 0.5 * Consts.sqrt3 * bc;
479 ydelta = y * Blob.delta;
480 dc = (REAL) fabs(ydelta - ay);
500 if (m >= lmf && m <= lml && b <= lbl)
502 if (n < Blob.support)
505 list[*numb] = h1 * Blob.M + m1;
506 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
508 weight[*numb] = Blob.ltable[nint];
509 *snorm += weight[*numb] * weight[*numb];
530 if (m >= lmf && m <= lml && b >= lbf)
532 if (n < Blob.support)
535 list[*numb] = h1 * Blob.M + m1;
536 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
538 weight[*numb] = Blob.ltable[nint];
540 *snorm += weight[*numb] * weight[*numb];
568 "\n bwray np = %5i nr = %5i numb = %5i snorm = %10.4f",
569 np, nr, *numb, *snorm);
572 if ((*numb > 0) && (trace > 9))
574 for (int i = 0; i < *numb; i++)
577 fprintf(output, "\n %5i %10.4f", list[i],
578 weight[i]); //bug 139