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/bray.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
26 void Blob_class::bray(INTEGER np, INTEGER nr, INTEGER* list, REAL* weight, INTEGER* numb, REAL* snorm)
30 INTEGER ac, bc, bc1, ac1;
43 REAL ax, ay, axid, ayid;
45 REAL theta, sinth, costh;
46 REAL prfx, prlx, prfy, prly;
49 REAL ddeltain, ddeltaout;
51 REAL ydelta; //xdelta not used
61 // bug 168 - swr - 9/24/05
65 "\n ***** bray is not usable with divergent geometry\n");
73 hex_voronoi = Consts.sqrt3 / 2.0 * Blob.delta * Blob.delta;
75 multf = GeoPar.pinc / (2.0 * Blob.delta); // pinc/2
76 pinc2 = GeoPar.pinc / 2.0;
78 laf = (INTEGER) (-0.5 * ((int) (Blob.M / 2) + (int) (Blob.H / 2)));
79 lal = (INTEGER) (0.5 * ((int) ((Blob.M - 1) / 2) + (int) ((Blob.H - 1) / 2)));
81 lbf = -1 * (INTEGER) (Blob.H / 2);
82 lbl = (INTEGER) ((Blob.H - 1) / 2);
84 lmf = -1 * (INTEGER) (Blob.M / 2);
85 lml = (INTEGER) ((Blob.M - 1) / 2);
87 prfx = -0.25 * (Blob.M - 1.0);
90 prfy = -0.25 * Consts.sqrt3 * (Blob.H - 1.0);
93 H2 = (INTEGER) (Blob.H / 2);
94 M2 = (INTEGER) (Blob.M / 2);
95 posit(np, nr, &ax, &ay, &mx, &my);
96 axid = ax / Blob.delta;
97 ayid = ay / Blob.delta;
104 theta = 2 * Consts.pi - theta; //theta is not from getang, but different from ray to ray, by wei 1/2005
106 //added 7/19, think of theta as always betw 0 and pi
107 if (theta > Consts.pi)
109 theta = theta - Consts.pi;
111 costh = -costh; //never used anyway
116 if (fabs(my) > Consts.zero)
117 { // non-horizontal lines
120 ord = (axid - islope * ayid);
121 coef = 0.5 * (Consts.sqrt3 * islope - 1);
123 if (theta < Consts.pi / 6 || theta > 2.0 * Consts.pi / 3.0)
126 cte = multf / my; // for the "side rays"
129 fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
131 mcov = (INTEGER) ceil(multf / fcosang);
132 ddeltain = Blob.delta * fcosang;
133 ddeltaout = Blob.delta * sinth;
135 if (theta < Consts.pi / 6)
146 // parameters of the first grid point
150 yup = (x - ord + ctes) * slope;
151 ylo = (x - ord - ctes) * slope;
159 x = ord - ctes + islope * yup;
161 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
166 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
170 x = ord + ctes + islope * ylo;
171 ac = (INTEGER) floor(x - ylo * Consts.isqrt3);
180 yup = (x - ord + ctes) * slope;
181 ylo = (x - ord - ctes) * slope;
189 x = ord + ctes + islope * yup;
191 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
196 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
200 x = ord - ctes + islope * yup;
201 ac = (INTEGER) ceil(x - yup * Consts.isqrt3);
207 //done with af and al
209 bc = (INTEGER) (floor(icoef * (a - ord)));
211 bf = MAX0(bc - mcov, lbf);
212 bl = MIN0(bc + mcov, lbl);
215 y = 0.5 * Consts.sqrt3 * bc;
216 dc = Blob.delta * (REAL) fabs(-my * (axid - x) + mx * (ayid - y));
227 //starting the loops--------------------------------
236 if (m >= lmf && m <= lml && b <= lbl)
241 list[*numb] = h1 * Blob.M + m1;
264 if (m >= lmf && m <= lml && b >= lbf)
269 list[*numb] = h1 * Blob.M + m1;
283 dc = dc + (REAL) sign * ddeltaout;
286 if (dc > ddeltain || dc < 0.0)
289 dc = (REAL) sign * (fabs(dc) - ddeltain);
292 bf = MAX0(bc - mcov, lbf);
293 bl = MIN0(bc + mcov, lbl);
309 fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
310 mcov = (INTEGER) ceil(multf / sinth);
312 // constants independent of the first grid point
315 ddeltain = Blob.delta * sinth;
316 ddeltaout = Blob.delta * fcosang;
318 if (theta < Consts.pi / 3)
327 //----------------------------------------------
328 // parameters of the first grid point
332 ac = (INTEGER) (floor(coef * b + ord));
334 af = MAX0(ac - mcov, laf);
335 al = MIN0(ac + mcov, lal);
338 y = 0.5 * Consts.sqrt3 * b;
339 dc = Blob.delta * (REAL) fabs(-my * (x - axid) + mx * (y - ayid));
359 if (m >= lmf && m <= lml)
364 list[*numb] = h1 * Blob.M + m1;
385 if (m >= lmf && m <= lml)
390 list[*numb] = h1 * Blob.M + m1;
402 dc = dc + (REAL) sign * ddeltaout;
405 if (dc > ddeltain || dc < 0.0)
409 dc = (REAL) sign * (fabs(dc) - ddeltain);
412 af = MAX0(ac - mcov, laf);
413 al = MIN0(ac + mcov, lal);
426 } // if non- horizontal line.
431 fcosang = cos(Consts.pi / 6.0);
433 mcov = (INTEGER) ceil(multf / fcosang);
435 // constants independent of the first grid point
436 ddeltain = Blob.delta * fcosang;
438 // parameters of the first grid point
443 ylo = yup - 2 * multf;
445 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
452 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
455 //done with af and al
458 bc = (INTEGER) (floor(2 * ayid * Consts.isqrt3));
459 bf = MAX0(bc - mcov, lbf);
460 bl = MIN0(bc + mcov, lbl);
462 y = 0.5 * Consts.sqrt3 * bc;
463 ydelta = y * Blob.delta;
464 dc = (REAL) fabs(ydelta - ay);
467 //nc1 = (INTEGER)(samp_pinc * dc1);
477 // starting the loops
486 if (m >= lmf && m <= lml && b <= lbl)
491 list[*numb] = h1 * Blob.M + m1;
513 if (m >= lmf && m <= lml && b >= lbf)
518 list[*numb] = h1 * Blob.M + m1;
545 for (i = 0; i < *numb; i++)
547 weight[i] = hex_voronoi;
551 *snorm = hex_voronoi * hex_voronoi * (*numb);
555 "\n bray np = %5i nr = %5i numb = %5i snorm = %10.4f",
556 np, nr, *numb, *snorm);
559 if ((*numb > 0) && (trace > 9))
561 for (int i = 0; i < *numb; i++)
564 fprintf(output, "\n %5i", list[i]); //bug 139