examples/b2/src/alp1.cpp

   1 /*
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                           *
   4  *                              S N A R K   1 4                              *
   5  *                                                                           *
   6  *                     A PICTURE RECONSTRUCTION PROGRAM                      *
   7  *                                                                           *
   8  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   9
  10  alp1.cpp,v 1.1 2009/06/01 23:26:59 jklukowska Exp
  11
  12 */
  13 //     Example 2  user-written reconstruction function
  14
  15
  16 #include <stdio.h>
  17 #include <math.h>
  18
  19 #include "blkdta.h"
  20 #include "pseudo.h"
  21 #include "projfile.h"
  22 #include "uiod.h"
  23 #include "consts.h"
  24
  25 #include "alp1.h"
  26
  27 void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area);
  28 void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight);
  29
  30 BOOLEAN alp1_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  31 {
  32
  33   static REAL* oldrec;
  34   static INTEGER curprj;
  35   static REAL* proj;
  36   INTEGER i;
  37   static REAL newdif;
  38   static REAL olddif;
  39   static REAL c;
  40
  41   if(iter == 1) {
  42
  43     oldrec = new REAL[GeoPar.area];
  44
  45     proj = new REAL[GeoPar.nrays];
  46
  47     for(i = 0; i < GeoPar.area; i++) {
  48       oldrec[i] = 0;
  49     }
  50
  51     // The previous art iteration is in OLDREC
  52     itrun(recon, oldrec, proj, list, weight);
  53
  54     for(i = 0; i < GeoPar.area; i++) {
  55       oldrec[i] = recon[i];
  56     }
  57
  58     curprj = 1;
  59     return FALSE;
  60   }
  61   else if(iter == 2) {
  62     itrun(recon, oldrec, proj, list, weight);
  63     newdif = 0.0;
  64     curprj = 2;
  65
  66     for(i = 0; i < GeoPar.area; i++) {
  67       newdif += SQR(recon[i] - oldrec[i]);
  68       oldrec[i] = recon[i];
  69     }
  70
  71     newdif = (REAL) sqrt(newdif);
  72     return FALSE;
  73   }
  74   else if((iter/2)*2 != iter) {
  75     curprj += 1;
  76     itrun(recon, oldrec, proj, list, weight);
  77     olddif = newdif;
  78     newdif = 0.0;
  79
  80     for(i = 0; i < GeoPar.area; i++) {
  81       newdif += SQR(recon[i] - oldrec[i]);
  82     }
  83
  84     newdif = (REAL) sqrt(newdif);
  85     if(olddif <= Consts.zero * newdif) {
  86       fprintf(output, "\nCONVERGENCE DIFFERENCES %20.12e, %20.12e", olddif, newdif);
  87       return FALSE;
  88     }
  89     c = newdif/olddif;
  90     fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
  91
  92     return FALSE;
  93   }
  94
  95   extrap(recon, oldrec, c, GeoPar.area);
  96
  97   return FALSE;
  98 }
  99
 100 // NEWREC will contain extrapolated values
 101 // OLDREC will contain the old values of NEWREC
 102
 103 void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area)
 104 {
 105   REAL d, save;
 106   INTEGER i;
 107
 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);
 111
 112     for(i = 0; i < area; i++) {
 113       save = newrec[i];
 114       newrec[i] += d * (newrec[i] - oldrec[i]);
 115       oldrec[i] = save;
 116     }
 117     return;
 118   }
 119   fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
 120
 121   fprintf(output, "\n EXTRAPOLATION PARAMETER TOO CLOSE TO 1 ");
 122
 123
 124   for(i = 0; i < area; i++) {
 125     oldrec[i] = newrec[i];
 126   }
 127   return;
 128 }
 129
 130 void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight)
 131 {
 132   INTEGER i, nr, np, numb, j;
 133   REAL snorm, raysum;
 134
 135   for(i = 0; i < GeoPar.area; i++) {
 136     recon[i] = oldrec[i];
 137   }
 138
 139   for(np = 0; np < GeoPar.prjnum; np++) {
 140
 141     ProjFile.ReadProj(np, proj, GeoPar.nrays);
 142
 143     for(nr = 0; nr < GeoPar.nrays; nr++) {
 144       raysum = pseudo(recon, np, nr, list, weight, &numb, &snorm, TRUE, FALSE);
 145
 146       for(i = 0; i < numb; i++) {
 147         j = list[i];
 148         recon[j] += (proj[nr]-raysum)* weight[i] * (REAL) 0.025 / snorm;
 149       }
 150     }
 151   }
 152
 153   return;
 154 }