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/bdhk_misl.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
38 REAL bdhk_class::AVER(REAL SUM, REAL NUM)
43 REAL bdhk_class::AVER3(REAL A, REAL B, REAL C)
45 REAL temp = (A + B + C) / 3.0;
49 REAL bdhk_class::diffAbs(REAL A, REAL B)
59 REAL bdhk_class::neigSum(REAL* recon, REAL* d)
62 REAL Pi1, Pi2, Pi3, Pi4, Dif1, Dif2, Dif3, Dif4;
63 for (int i = 0; i < GeoPar.area; i++)
65 if (i == 0) //upper left corner
67 Pi1 = AVER3(recon[1], recon[GeoPar.nelem], recon[GeoPar.nelem + 1]);
68 Dif1 = diffAbs(recon[i], Pi1);
70 d[i] = Pi1 - recon[i];
72 else if (i == (GeoPar.nelem - 1)) //upper right corner
74 Pi1 = AVER3(recon[GeoPar.nelem - 2], recon[2 * GeoPar.nelem - 2],
75 recon[2 * GeoPar.nelem - 1]);
76 Dif1 = diffAbs(recon[i], Pi1);
78 d[i] = Pi1 - recon[i];
80 else if (i == (GeoPar.area - GeoPar.nelem)) //lower left corner
82 Pi1 = AVER3(recon[GeoPar.area - 2 * GeoPar.nelem],
83 recon[GeoPar.area - 2 * GeoPar.nelem + 1],
84 recon[GeoPar.area - GeoPar.nelem + 1]);
85 Dif1 = diffAbs(recon[i], Pi1);
87 d[i] = Pi1 - recon[i];
89 else if (i == (GeoPar.area - 1)) //lower right corner
91 Pi1 = AVER3(recon[GeoPar.area - 2],
92 recon[GeoPar.area - GeoPar.nelem - 1],
93 recon[GeoPar.area - GeoPar.nelem - 2]);
94 Dif1 = diffAbs(recon[i], Pi1);
96 d[i] = Pi1 - recon[i];
99 else if (i < GeoPar.nelem) //first row, not corners
103 Pi1 = AVER3(recon[i - 1], recon[i - 1 + GeoPar.nelem],
104 recon[i + GeoPar.nelem]);
105 Pi2 = AVER3(recon[i + 1], recon[i + 1 + GeoPar.nelem],
106 recon[i + GeoPar.nelem]);
108 Dif1 = diffAbs(recon[i], Pi1);
109 Dif2 = diffAbs(recon[i], Pi2);
114 d[i] = Pi1 - recon[i];
120 d[i] = Pi2 - recon[i];
124 else if (i % GeoPar.nelem == 0) //first column, not corners
128 Pi1 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
129 recon[i - GeoPar.nelem + 1]);
130 Pi2 = AVER3(recon[i + 1], recon[i + GeoPar.nelem],
131 recon[i + GeoPar.nelem + 1]);
133 Dif1 = diffAbs(recon[i], Pi1);
134 Dif2 = diffAbs(recon[i], Pi2);
139 d[i] = Pi1 - recon[i];
144 d[i] = Pi2 - recon[i];
149 else if (i % GeoPar.nelem == GeoPar.nelem - 1) //last column, not corners
153 Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
154 recon[i - GeoPar.nelem]);
155 Pi2 = AVER3(recon[i - 1], recon[i + GeoPar.nelem],
156 recon[i + GeoPar.nelem - 1]);
158 Dif1 = diffAbs(recon[i], Pi1);
159 Dif2 = diffAbs(recon[i], Pi2);
164 d[i] = Pi1 - recon[i];
169 d[i] = Pi2 - recon[i];
173 else if (i > GeoPar.area - GeoPar.nelem) //last row, not corners
177 Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
178 recon[i - GeoPar.nelem]);
179 Pi2 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
180 recon[i - GeoPar.nelem + 1]);
182 Dif1 = diffAbs(recon[i], Pi1);
183 Dif2 = diffAbs(recon[i], Pi2);
188 d[i] = Pi1 - recon[i];
193 d[i] = Pi2 - recon[i];
197 else //not in the bounderies but in the middle
199 //four possibilities:
201 Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
202 recon[i - GeoPar.nelem]);
203 Pi2 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
204 recon[i - GeoPar.nelem + 1]);
205 Pi3 = AVER3(recon[i - 1], recon[i - 1 + GeoPar.nelem],
206 recon[i + GeoPar.nelem]);
207 Pi4 = AVER3(recon[i + 1], recon[i + GeoPar.nelem],
208 recon[i + GeoPar.nelem + 1]);
210 Dif1 = diffAbs(recon[i], Pi1);
211 Dif2 = diffAbs(recon[i], Pi2);
212 Dif3 = diffAbs(recon[i], Pi3);
213 Dif4 = diffAbs(recon[i], Pi4);
215 if (Dif1 <= Dif2 && Dif1 <= Dif3 && Dif1 <= Dif4)
218 d[i] = Pi1 - recon[i];
220 else if (Dif2 <= Dif1 && Dif2 <= Dif3 && Dif2 <= Dif4)
223 d[i] = Pi2 - recon[i];
225 else if (Dif3 <= Dif1 && Dif3 <= Dif2 && Dif3 <= Dif4)
228 d[i] = Pi3 - recon[i];
230 else if (Dif4 <= Dif1 && Dif4 <= Dif2 && Dif4 <= Dif3)
233 d[i] = Pi4 - recon[i];
238 REAL norm_d = calcNorm(d);
242 //normalizing the d vector:
243 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~
244 for (int i1 = 0; i1 < GeoPar.area; i1++)
251 REAL bdhk_class::calculateResidual(REAL* recon, INTEGER* list, REAL* weight,
252 struct NPNR* arrayRays, INTEGER kount)
254 INTEGER np, nr, i, j, numb;
255 REAL snorm, temp1, b, res = 0;
256 for (i = 0; i < kount; i++)
258 np = arrayRays[i].NP;
259 nr = arrayRays[i].NR;
260 wray(np, nr, list, weight, &numb, &snorm);
261 b = Anglst.prdta(np, nr);
265 for (j = 0; j < numb; j++)
266 temp1 += weight[j] * recon[list[j]];
267 res += (pow((b - temp1), 2)) / snorm;
274 BOOLEAN bdhk_class::calculateGradientTV(REAL* grad, REAL* recon)
276 //initialiaze the Gradient vector to all zero
277 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
280 //looping through case 1
281 for (int k = 0; k < (GeoPar.area - GeoPar.nelem); k++) //not including the last row
283 if (((k % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not in the last column
286 REAL v1 = recon[k + 1] - recon[k];
287 REAL v2 = recon[k + GeoPar.nelem] - recon[k];
288 REAL deno = sqrt(v1 * v1 + v2 * v2);
289 if (deno >= Consts.zero) //continue calculating the subgradients
290 //there are three subgradients to consider. They are placed in k, k+1, k+nelem
291 // (1) 2x_k - x_k+1 - x_k+nelem / deno
292 grad[k] += (2 * recon[k] - recon[k + 1]
293 - recon[k + GeoPar.nelem]) / deno;
299 //looping through case 2
300 for (int k = 0; k < (GeoPar.area - GeoPar.nelem); k++) //not including the last row
302 if ((k % GeoPar.nelem) != 0) //not in the first column
305 REAL v1 = recon[k] - recon[k - 1];
306 REAL v2 = recon[k + GeoPar.nelem - 1] - recon[k - 1];
307 REAL deno = sqrt(v1 * v1 + v2 * v2);
308 if (deno >= Consts.zero) //continue calculating the subgradients
309 grad[k] += (recon[k] - recon[k - 1]) / deno;
315 //looping through case 3
316 for (int k = GeoPar.nelem; k < (GeoPar.area); k++) //not including the first row
318 if (((k % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not in the last column
321 REAL v1 = recon[k + 1 - GeoPar.nelem] - recon[k - GeoPar.nelem];
322 REAL v2 = recon[k] - recon[k - GeoPar.nelem];
323 REAL deno = sqrt(v1 * v1 + v2 * v2);
324 if (deno >= Consts.zero) //continue calculating the subgradients
325 grad[k] += (recon[k] - recon[k - GeoPar.nelem]) / deno;
331 // BOOLEAN flag=FALSE;
332 for (int i = 0; i < GeoPar.area; i++)
340 } //end calculateGradient
342 void bdhk_class::orderRays(struct NPNR* arrayRays, INTEGER kount)
345 for (i = 0; i < kount; i++)
355 void bdhk_class::MyART(REAL* recon, INTEGER* list, REAL* weight,
356 struct NPNR* arrayRays, INTEGER kount)
358 INTEGER np, nr, i, j, numb;
359 REAL snorm, temp, b, div, dif;
361 for (i = 0; i < kount; i++)
363 np = arrayRays[i].NP;
364 nr = arrayRays[i].NR;
365 wray(np, nr, list, weight, &numb, &snorm);
368 if (snorm <= Consts.zero)
370 cout << endl << "*** error ---> norm square is less than "
371 << Consts.zero << " ***";
372 cout << endl << " *** ART CALCULATION ABORTED ***" << endl;
376 for (j = 0; j < numb; j++)
377 temp += weight[j] * recon[list[j]];
378 b = Anglst.prdta(np, nr);
381 for (j = 0; j < numb; j++)
382 recon[list[j]] += div * weight[j];
389 REAL bdhk_class::calcHalfNormSquare(REAL* recon)
391 REAL sum = squareNorm(recon) / 2;
395 REAL bdhk_class::calcTV(REAL* recon)
399 for (int x = 0; x < (GeoPar.area - GeoPar.nelem - 1); x++)
401 if (((x % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not last column
403 pow((recon[x + 1] - recon[x]), 2)
404 + pow((recon[x + GeoPar.nelem] - recon[x]), 2));
411 REAL bdhk_class::squareNorm(REAL* aVector)
414 for (int i = 0; i < GeoPar.area; i++)
415 sum += aVector[i] * aVector[i];
419 REAL bdhk_class::calcNorm(REAL* aVector)
421 return sqrt(squareNorm(aVector));