2 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
6 * A PICTURE RECONSTRUCTION PROGRAM *
8 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
10 alp1.cpp,v 1.1 2009/06/01 23:26:59 jklukowska Exp
13 // Example 2 user-written reconstruction function
27 void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area);
28 void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight);
30 BOOLEAN alp1_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
34 static INTEGER curprj;
43 oldrec = new REAL[GeoPar.area];
45 proj = new REAL[GeoPar.nrays];
47 for(i = 0; i < GeoPar.area; i++) {
51 // The previous art iteration is in OLDREC
52 itrun(recon, oldrec, proj, list, weight);
54 for(i = 0; i < GeoPar.area; i++) {
62 itrun(recon, oldrec, proj, list, weight);
66 for(i = 0; i < GeoPar.area; i++) {
67 newdif += SQR(recon[i] - oldrec[i]);
71 newdif = (REAL) sqrt(newdif);
74 else if((iter/2)*2 != iter) {
76 itrun(recon, oldrec, proj, list, weight);
80 for(i = 0; i < GeoPar.area; i++) {
81 newdif += SQR(recon[i] - oldrec[i]);
84 newdif = (REAL) sqrt(newdif);
85 if(olddif <= Consts.zero * newdif) {
86 fprintf(output, "\nCONVERGENCE DIFFERENCES %20.12e, %20.12e", olddif, newdif);
90 fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
95 extrap(recon, oldrec, c, GeoPar.area);
100 // NEWREC will contain extrapolated values
101 // OLDREC will contain the old values of NEWREC
103 void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area)
108 if(fabs(1.0 - c) > Consts.zero) {
109 fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
110 d = c / ((REAL) 1.0 - c);
112 for(i = 0; i < area; i++) {
114 newrec[i] += d * (newrec[i] - oldrec[i]);
119 fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
121 fprintf(output, "\n EXTRAPOLATION PARAMETER TOO CLOSE TO 1 ");
124 for(i = 0; i < area; i++) {
125 oldrec[i] = newrec[i];
130 void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight)
132 INTEGER i, nr, np, numb, j;
135 for(i = 0; i < GeoPar.area; i++) {
136 recon[i] = oldrec[i];
139 for(np = 0; np < GeoPar.prjnum; np++) {
141 ProjFile.ReadProj(np, proj, GeoPar.nrays);
143 for(nr = 0; nr < GeoPar.nrays; nr++) {
144 raysum = pseudo(recon, np, nr, list, weight, &numb, &snorm, TRUE, FALSE);
146 for(i = 0; i < numb; i++) {
148 recon[j] += (proj[nr]-raysum)* weight[i] * (REAL) 0.025 / snorm;