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_deset.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
14 COPT = 2 MODIFIED SIRT
34 void quad_class::deset(REAL* d, INTEGER* list, REAL* weight, BOOLEAN* alg)
57 for (i = 0; i < area; i++)
62 for (np = 0; np < GeoPar.prjnum; np++)
64 for (nr = 0; nr < GeoPar.nrays; nr++)
69 wray(np, nr, list, weight, &numb, &snorm);
71 Blob.bwray(np, nr, list, weight, &numb, &snorm);
76 ray(np, nr, list, weight, &numb, &snorm);
78 Blob.bray(np, nr, list, weight, &numb, &snorm);
89 for (k = 0; k < numb; k++)
98 truval = Anglst.prdta(np, nr);
99 error_1 = uerror(np, nr, truval, alg);
103 truval = Anglst.prdta(np, nr);
104 error_1 = errfac(truval, 0.0, 0.0);
109 for (k = 0; k < numb; k++)
112 d[j] += totwei * weight[k];
123 truval = Anglst.prdta(np, nr);
124 error_1 = uerror(np, nr, truval, alg);
128 truval = Anglst.prdta(np, nr);
129 error_1 = errfac(truval, 0.0, 0.0);
133 for (k = 0; k < numb; k++)
136 d[j] += (weight[k] * weight[k]) / error_1;
148 for (i = 0; i < area; i++)
150 if (fabs(d[i]) <= Consts.zero)
153 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
156 fprintf(output, "\n\n ***a picture element is missed by all rays, can not compute matrix D, program stops***"); //hstau
161 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
163 d[i] = (REAL) 1.0 / d[i];