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_mtamxp.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
12 X AND Y MUST BE DISTINCT
31 void quad_class::mtamxp(REAL* x, REAL* y, INTEGER* list, REAL* weight)
52 hex_voronoi = Consts.sqrt3 / (REAL) 2.0 * Blob.delta * Blob.delta;
67 //for(i = 0; i < area; i++) {
68 // fprintf(output,"\nmtamxpx=%f", x[i]);
74 for (i = 0; i < area; i++)
79 // FOR BETTER PERFORMANCE LINE AND STRIP MODES ARE SEPERATED
86 for (np = 0; np < GeoPar.prjnum; np++)
88 for (nr = 0; nr < GeoPar.nrays; nr++)
93 wray(np, nr, list, weight, &numb, &snorm);
97 Blob.bwray(np, nr, list, weight, &numb, &snorm);
105 for (nb = 0; nb < numb; nb++)
108 val += x[k] * weight[nb];
113 truval = Anglst.prdta(np, nr);
116 // VAL = SA*A*(M*X - P)
122 error_1 = uerror(np, nr, truval, &alg);
126 error_1 = errfac(truval, 0.0, 0.0);
132 // Y = SA*M#*A*(M*X - P)
134 for (nb = 0; nb < numb; nb++)
137 y[k] += val * weight[nb];
152 Fourie.rinc = GeoPar.pinc; //think rinc and these pix* for blob case. hstau
155 pix2 = GeoPar.pixsiz * GeoPar.pixsiz;
157 pix4 = pix3; // bug in Fortran?. NO hstau
166 for (np = 0; np < GeoPar.prjnum; np++)
168 Anglst.getang(np, &theta, &sinth, &costh);
169 if ((NoisePar.quanin > 0) && GeoPar.vri)
170 Fourie.rinc = GeoPar.pinc
171 * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
172 rinc2 = Fourie.rinc * Fourie.rinc;
174 for (nr = 0; nr < GeoPar.nrays; nr++)
178 ray(np, nr, list, weight, &numb, &snorm);
182 Blob.bray(np, nr, list, weight, &numb, &snorm);
189 // VAL = M*X (NOTE THAT M IS A 0-1 MATRIX)
191 for (nb = 0; nb < numb; nb++)
201 truval = Anglst.prdta(np, nr);
205 // VAL = SA*M#*A*(M*X - P)
211 error_1 = uerror(np, nr, truval, &alg);
215 error_1 = errfac(truval, Fourie.rinc, rinc2);
218 val *= pix4 / error_1;
220 // Y = SA*M#*A*(M*X - P)
222 for (nb = 0; nb < numb; nb++)