Added snark14m distribution examples

[snark14.git] / examples / b2 / src / alp1.cpp

diff --git a/examples/b2/src/alp1.cpp b/examples/b2/src/alp1.cpp
new file mode 100644 (file)
index 0000000..c4835a3
--- /dev/null
+++ b/examples/b2/src/alp1.cpp
@@ -0,0 +1,154 @@
+/* 
+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ *                                                                           *
+ *                              S N A R K   1 4                              *
+ *                                                                           *
+ *                     A PICTURE RECONSTRUCTION PROGRAM                      *
+ *                                                                           *
+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+ alp1.cpp,v 1.1 2009/06/01 23:26:59 jklukowska Exp
+
+*/
+//     Example 2  user-written reconstruction function
+
+
+#include <stdio.h>
+#include <math.h>
+
+#include "blkdta.h"
+#include "pseudo.h"
+#include "projfile.h"
+#include "uiod.h"
+#include "consts.h"
+
+#include "alp1.h"
+
+void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area);
+void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight);
+
+BOOLEAN alp1_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
+{
+
+  static REAL* oldrec;
+  static INTEGER curprj;
+  static REAL* proj;
+  INTEGER i;
+  static REAL newdif; 
+  static REAL olddif;
+  static REAL c;
+
+  if(iter == 1) {
+
+    oldrec = new REAL[GeoPar.area];
+
+    proj = new REAL[GeoPar.nrays];
+
+    for(i = 0; i < GeoPar.area; i++) {
+      oldrec[i] = 0;
+    }
+
+    // The previous art iteration is in OLDREC
+    itrun(recon, oldrec, proj, list, weight);
+
+    for(i = 0; i < GeoPar.area; i++) {
+      oldrec[i] = recon[i];
+    }
+
+    curprj = 1;
+    return FALSE;
+  }
+  else if(iter == 2) {
+    itrun(recon, oldrec, proj, list, weight);
+    newdif = 0.0;
+    curprj = 2;
+
+    for(i = 0; i < GeoPar.area; i++) {
+      newdif += SQR(recon[i] - oldrec[i]);
+      oldrec[i] = recon[i];
+    }
+
+    newdif = (REAL) sqrt(newdif);
+    return FALSE;
+  }
+  else if((iter/2)*2 != iter) {
+    curprj += 1;
+    itrun(recon, oldrec, proj, list, weight);
+    olddif = newdif;
+    newdif = 0.0;
+
+    for(i = 0; i < GeoPar.area; i++) {
+      newdif += SQR(recon[i] - oldrec[i]);
+    }
+
+    newdif = (REAL) sqrt(newdif);
+    if(olddif <= Consts.zero * newdif) {
+      fprintf(output, "\nCONVERGENCE DIFFERENCES %20.12e, %20.12e", olddif, newdif);
+      return FALSE;
+    }
+    c = newdif/olddif;
+    fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
+
+    return FALSE;
+  }
+
+  extrap(recon, oldrec, c, GeoPar.area);
+
+  return FALSE;
+}
+
+// NEWREC will contain extrapolated values
+// OLDREC will contain the old values of NEWREC
+
+void extrap(REAL* newrec, REAL* oldrec, REAL c, INTEGER area)
+{
+  REAL d, save;
+  INTEGER i;
+
+  if(fabs(1.0 - c) > Consts.zero) {
+    fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
+    d = c / ((REAL) 1.0 - c);
+
+    for(i = 0; i < area; i++) {
+      save = newrec[i];
+      newrec[i] += d * (newrec[i] - oldrec[i]);
+      oldrec[i] = save;
+    }
+    return;
+  }
+  fprintf(output, "\n EXTRAPOLATION PARAMETER IN EXTRAP IS %20.12e", c);
+
+  fprintf(output, "\n EXTRAPOLATION PARAMETER TOO CLOSE TO 1 ");
+
+
+  for(i = 0; i < area; i++) {
+    oldrec[i] = newrec[i];
+  }
+  return;
+}
+
+void itrun(REAL* recon, REAL* oldrec, REAL* proj, INTEGER* list, REAL* weight)
+{
+  INTEGER i, nr, np, numb, j;
+  REAL snorm, raysum;
+
+  for(i = 0; i < GeoPar.area; i++) {
+    recon[i] = oldrec[i];
+  }
+
+  for(np = 0; np < GeoPar.prjnum; np++) {
+
+    ProjFile.ReadProj(np, proj, GeoPar.nrays);
+
+    for(nr = 0; nr < GeoPar.nrays; nr++) {
+      raysum = pseudo(recon, np, nr, list, weight, &numb, &snorm, TRUE, FALSE);
+
+      for(i = 0; i < numb; i++) {
+        j = list[i];
+        recon[j] += (proj[nr]-raysum)* weight[i] * (REAL) 0.025 / snorm;
+      }
+    }
+  }
+
+  return;
+}