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_costfn.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 THIS SUBROUTINE COMPUTES THE VALUE OF K(X) NEEDED BY THE QUADRATIC
11 OPTIMIZATION ALGORITHM.
30 void quad_class::costfn(REAL* recon, REAL* r, REAL* g, REAL* cg, REAL* w,
31 INTEGER* list, REAL* weight, REAL* f)
34 static REAL ppcost; // made static to avoid garbage value for `iterx > 1'. Lajos, Feb 10, 2005
68 fprintf(output, "\n\n cost function can not be computed");
75 Fourie.rinc = GeoPar.pinc;
76 for (np = 0; np < GeoPar.prjnum; np++)
78 Anglst.getang(np, &theta, &sinth, &costh);
80 Fourie.rinc = GeoPar.pinc
81 * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
82 rinc2 = Fourie.rinc * Fourie.rinc;
83 for (nr = 0; nr < GeoPar.nrays; nr++)
85 pr = Anglst.prdta(np, nr);
90 error_1 = uerror(np, nr, pr, &alg);
94 error_1 = errfac(pr, Fourie.rinc, rinc2);
97 ppcost += pr * pr / error_1;
108 if (!((sb < Consts.zero) || (l <= 0)))
120 if ((sb + sc) > Consts.zero)
122 for (i = 0; i < area; i++)
126 n[i] = recon[i] - xbar;
130 if ((sb > Consts.zero) && (l > 0))
134 adsmos(n, l, bw1, bw2, bw3, bbcon);
138 badsmos(n, l, bw1, bw2, bbcon);
142 for (i = 0; i < area; i++)
144 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
146 xxcost += m[i] * m[i];
147 if (sb > Consts.zero)
148 bxcost += m[i] * n[i];
155 Fourie.rinc = GeoPar.pinc;
156 for (np = 0; np < GeoPar.prjnum; np++)
158 Anglst.getang(np, &theta, &sinth, &costh);
160 Fourie.rinc = GeoPar.pinc * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
161 rinc2 = Fourie.rinc * Fourie.rinc;
162 for (nr = 0; nr < GeoPar.nrays; nr++)
164 pr = Anglst.prdta(np, nr);
169 ray(np, nr, list, weight, &numb, &snorm);
173 Blob.bray(np, nr, list, weight, &numb, &snorm);
180 wray(np, nr, list, weight, &numb, &snorm);
184 Blob.bwray(np, nr, list, weight, &numb, &snorm);
190 for (nb = 0; nb < numb; nb++)
193 val += recon[k] * weight[nb];
200 error_1 = uerror(np, nr, pr, &alg);
204 error_1 = errfac(pr, Fourie.rinc, rinc2);
207 valder = val / error_1;
208 xmcost += val * valder;
209 pmcost += pr * valder;
214 cost = sa * (ppcost - (REAL) 2.0 * pmcost + xmcost) + sb * bxcost + sc * xxcost;
215 fprintf(output, "\n\n value of cost function %15.6e", cost);