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/quad_badsmos.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 THIS ROUTINE PERFORMS THE FOUR DIFFERENT SMOOTHING PROCEDURES
11 USED IN THE QUADRATIC OPTIMIZATION ALGORITHM FOR BLOB RECONSTRUCTIONS.//changed . hstau
23 //pragma message( "to fix!" )
24 /// fix indexes 1 -> 0
26 void quad_class::badsmos(REAL* x, INTEGER m, REAL w1, REAL w2, INTEGER bcon)
30 static const INTEGER hsnar = CHAR2INT('s', 'n', 'a', 'r');
31 static const INTEGER hsajc = CHAR2INT('s', 'a', 'j', 'c');
32 static const INTEGER hkami = CHAR2INT('k', 'a', 'm', 'i');
35 //INTEGER n, // number of rows
36 INTEGER H; // number of rows
37 INTEGER M; // number of columns
38 INTEGER m1; // Blob.M-1
39 INTEGER m2; //Blob.M+1
40 INTEGER h1; //Blob.H-1
51 REAL* in = NULL; //wei 3/2005
52 //REAL* ipf; //wei 3/2005
69 // FOR THE KAMI METHOD THE ROUTINE IS EXECUTED TWICE
92 weit4 = (REAL) (w1 + w2 + 2.0 * w3);
93 weit6 = weit4 + w2 + w3;
94 //weit9 = (REAL) (weit6 + w2 + 2.0 * w3);
95 weit3 = (REAL) (w1 + w2 + w3);
96 weit5 = weit3 + w2 + w3;
97 weit7 = (REAL) (weit5 + 2.0 * w3);
100 // ALLOCATE WORKING SPACE
110 for (j = 0; j < m; j++)
113 // SMOOTH FIRST ROW OF PICTURE
115 in[Blob.pr] = x[Blob.pr];
116 in[m1 - Blob.pr] = x[m1 - Blob.pr];
122 x[0] = w1 * x[0] + w2 * x[2] + w3 * x[m2];
124 x[Blob.pr] /= weit3; //x[0] /= weit4;
129 x[1] = w1 * x[1] + w2 * x[3] + w3 * x[Blob.pr + m2]
132 x[Blob.pr] /= weit4; //x[0] /= weit4;
134 for (i = Blob.pr + 2; i < m1 - Blob.pr; i += 2)
137 x[i] = w1 * x[i] + w2 * (in[i - 2] + x[i + 2])
138 + w3 * (x[m1 + i] + x[m2 + i]);
141 x[i] /= weit5; //x[i] /= weit6;
145 x[m1] = w1 * x[m1] + w2 * in[m1 - 2] + w3 * x[m1 + m1];
146 //if(snar) x[n-1] /= weit4;
152 x[m1 - 1] = w1 * x[m1 - 1] + w2 * in[m1 - 3]
153 + w3 * x[m1 - 1 + m1] + w3 * x[m1 - 1 + m2];
165 REAL a2 = (x[0] + x[m2] + x[2 * M]) / 3.0;
166 x[0] = (REAL) (w1 * x[0] + w2 * (a1 + x[2])
167 + w3 * (x[0] + (x[0] + x[2]) / 2.0 + a2 + x[m2]));
171 REAL a1 = (x[1] + x[M]) / 2.0;
172 x[1] = (REAL) (w1 * x[1] + w2 * (a1 + x[3])
174 * (x[1] + (x[1] + x[3]) / 2.0 + x[M]
178 for (i = Blob.pr + 2; i < m1 - Blob.pr; i += 2)
181 x[i] = w1 * x[i] + w2 * (in[i - 2] + x[i + 2]) + w3 * (x[m1 + i] + x[m2 + i] + x[i] + (in[i - 2] + x[i + 2]) / 2.0);
185 // bug 191 - SAJC smoothing - swr - 12/17/05
189 REAL a2 = (x[m1] + x[m1 + m1] + x[2 * M + m1]) / 3.0;
190 x[m1] = (REAL) (w1 * x[m1] + w2 * (in[m1 - 2] + a1) + w3 * (x[m1] + (x[m1] + in[m1 - 2]) / 2.0 + x[m1 + m1] + a2));
194 REAL a1 = (x[m1 - 1] + x[m1 - 1 + m2]) / 2.0;
195 x[m1 - 1] = (REAL) (w1 * x[m1 - 1] + w2 * (in[m1 - 3] + a1) + w3 * ((in[m1 - 3] + x[m1 - 1]) / 2.0 + x[m1 - 1] + x[m1 + m1 - 1] + x[m1 + m1 + 1]));
200 for (k = 1; k < h1; k++)
205 kpr = (k - Blob.pr) % 2;
208 ln = l1 + m1 - 2 * kpr;
213 in[m1 - kpr] = x[ln];
215 x[l1] = w1 * x[l1] + w2 * x[l1 + 2];
219 x[l1] += w3 * (ip[1] + x[l1 + m2]);
223 x[l1] += w3 * (ip[0] + ip[2] + x[l1 + m2] + x[l1 + m1]);
227 // bug 191 - SAJC smoothing - swr - 12/17/05
232 a1 = (lt + ip[1] + in[0]) / 3.0;
234 a1 = (ip[1] + in[0]) / 2.0;
238 a3 = (in[0] + x[l1 + m2] + x[l1 + 2 * M]) / 3.0;
240 a3 = (in[0] + x[l1 + m2]) / 2.0;
241 x[l1] += w2 * a2 + w3 * (a1 + a3);
245 REAL a1 = (ip[0] + in[1] + x[l1 + m1]) / 3.0;
257 for (i = 2 + kpr; i < m1 - kpr; i += 2)
258 { //next point on the same row
261 x[li] = w1 * x[li] + w2 * (in[i - 2] + x[li + 2]) + w3 * (ip[i - 1] + ip[i + 1] + x[li + m1] + x[li + m2]);
266 x[ln] = w1 * x[ln] + w2 * in[m1 - 2 - kpr];
270 x[ln] += w3 * (ip[m1 - 1] + x[ln + m1]);
275 * (ip[m1 - 2] + x[ln + m1] + x[ln + m2] + ip[m1]);
288 // bug 191 - SAJC smoothing - swr - 12/17/05
293 a1 = (rt + in[m1] + ip[m1 - 1]) / 3.0;
295 a1 = (in[m1] + ip[m1 - 1]) / 2.0;
299 a3 = (in[m1] + x[ln + m1] + x[ln + 2 * M]) / 3.0;
301 a3 = (in[m1] + x[ln + m1]) / 2.0;
302 x[ln] += w2 * a2 + w3 * (a1 + a3);
306 REAL a1 = (in[m1 - 1] + ip[m1] + x[ln + m2]) / 3.0;
314 // IN LAST ROW CHANGE ROLE OF IN AND IP
316 l1 = h1 * M + Blob.pr;
317 ln = l1 + m1 - 2 * Blob.pr;
320 rt = ip[m1 - Blob.pr];
328 x[l1] = w1 * x[l1] + w2 * x[l1 + 2] + w3 * (in[1]);
332 x[l1] = w1 * x[l1] + w2 * x[l1 + 2] + w3 * (in[0] + in[2]);
343 for (i = 2 + Blob.pr; i < m1 - Blob.pr; i += 2)
348 x[li] = w1 * x[li] + w2 * (ip[i - 2] + x[li + 2])
349 + w3 * (in[i - 1] + in[i + 1]);
356 x[ln] = w1 * x[ln] + w2 * ip[m1 - 2] + w3 * in[m1 - 1];
360 x[ln] = w1 * x[ln] + w2 * ip[m1 - 3]
361 + w3 * (in[m1 - 2] + in[m1]);
374 // bug 191 - SAJC smoothing - swr - 12/17/05
378 REAL a2 = (x[l1] + lt + in[1]) / 3.0;
379 x[l1] = w1 * x[l1] + w2 * (x[l1 + 2] + a1)
381 * (a2 + in[1] + x[l1]
382 + (x[l1] + x[l1 + 2]) / 2.0);
386 REAL a1 = (x[l1] + in[0]) / 2.0;
387 x[l1] = w1 * x[l1] + w2 * (a1 + x[l1 + 2]) + w3 * (in[0] + in[2] + x[l1] + (x[l1] + x[l1 + 2]) / 2.0);
390 for (i = 2 + Blob.pr; i < m1 - Blob.pr; i += 2)
395 x[li] = w1 * x[li] + w2 * (ip[i - 2] + x[li + 2])
397 * (x[li] + 0.5 * ip[i - 2] + 0.5 * x[li + 2]
398 + in[i - 1] + in[i + 1]);
402 // bug 191 - SAJC smoothing - swr - 12/17/05
406 REAL a2 = (x[ln] + in[m1 - 1] + rt) / 3.0;
407 x[ln] = w1 * x[ln] + w2 * (ip[m1 - 2] + a1) + w3 * (a2 + in[m1 - 1] + (ip[m1 - 2] + x[ln]) / 2.0 + x[ln]);
411 REAL a1 = (x[ln] + in[m1]) / 2.0;
412 x[ln] = w1 * x[ln] + w2 * (ip[m1 - 3] + a1) + w3 * (in[m1 - 2] + in[m1] + (ip[m1 - 3] + x[ln]) / 2.0 + x[ln]);
417 if ((bcon != hkami) || (kami == 2))
420 // PREPARE FOR SECOND ROUND FOR KAMI
424 for (i = 0; i < M; i++)
433 for (i = 1; i < h1; i++)
443 // delete[] ipf; // bug 92 - Lajos - 03/02/2005 wei/2005
445 delete[] in; // wei 3/2005