r114: *** empty log message ***

[ctsim.git] / libctsim / phantom.cpp

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index 1508bea745912e3ffc07661c057f52c2096db4b1..ec5d0dcf6e0525cc38af1f1365091dbb81889e4b 100644 (file)
--- a/libctsim/phantom.cpp
+++ b/libctsim/phantom.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phantom.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**  $Id: phantom.cpp,v 1.3 2000/06/19 20:08:09 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -58,13 +58,13 @@ Phantom::create (const int phmid)
   switch (phmid) 
     {
     case O_PHM_HERMAN:
-      std_herman();
+      addStdHerman();
       break;
     case O_PHM_ROWLAND:
-      std_rowland();
+      addStdRowland();
       break;
     case O_PHM_BROWLAND:
-      std_rowland_bordered ();
+      addStdRowlandBordered ();
       break;
     case O_PHM_UNITPULSE:
       m_composition = P_UNIT_PULSE;
@@ -238,21 +238,18 @@ Phantom::show (void) const
  *   draw ()
  */
 
-#ifdef HAVE_SGP
 void 
 Phantom::draw (void) const
 {
+#ifdef HAVE_SGP
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
     sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
-}
 #endif
+}
 
 
 /* NAME
- *   std_rowland               Make head phantom of S.W. Rowland
- *
- * SYNOPSIS
- *   std_rowland ()
+ *   addStdRowland             Make head phantom of S.W. Rowland
  *
  * REFERENCES
  *   S. W. Rowland, "Computer Implementation of Image Reconstruction
@@ -261,7 +258,7 @@ Phantom::draw (void) const
  */
 
 void 
-Phantom::std_rowland (void)
+Phantom::addStdRowland (void)
 {
   addPElem("ellipse",  0.0000,  0.0000, 0.6900,  0.9200,   0.0,  1.00);
   addPElem("ellipse",  0.0000, -0.0184, 0.6624,  0.8740,   0.0, -0.98);
@@ -277,17 +274,14 @@ Phantom::std_rowland (void)
 }
 
 void 
-Phantom::std_rowland_bordered (void)
+Phantom::addStdRowlandBordered (void)
 {
-  std_rowland ();
+  addStdRowland ();
   addPElem ("ellipse", 0.000, 0.0000, 0.7500, 1.000, 0.0, 0.00);
 }
 
 /* NAME
- *   std_herman                        Standard head phantom of G. T. Herman
- *
- * SYNOPSIS
- *   std_herman ()
+ *   addStdHerman                      Standard head phantom of G. T. Herman
  *
  * REFERENCES
  *   G. T. Herman, "Image Reconstructions from Projections:  The Fundementals
@@ -295,7 +289,7 @@ Phantom::std_rowland_bordered (void)
  */
 
 void 
-Phantom::std_herman (void)
+Phantom::addStdHerman (void)
 {
   addPElem("ellipse",  0.000,  1.50,  0.375, 0.3000,  90.00, -0.003);
   addPElem("ellipse",  0.675, -0.75,  0.225, 0.1500, 140.00,  0.010);
@@ -315,6 +309,90 @@ Phantom::std_herman (void)
 }
 
 
+/* NAME
+ *    convertToImagefile               Make image array from Phantom
+ *
+ * SYNOPSIS
+ *    pic_to_imagefile (pic, im, nsample)
+ *    Phantom& pic             Phantom definitions
+ *    ImageFile  *im           Computed pixel array
+ *    int nsample              Number of samples along each axis for each pixel
+ *                             (total samples per pixel = nsample * nsample)
+ */
+
+void
+Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const
+{
+  convertToImagefile (im, in_nsample, trace, 0, im.nx());
+}
+
+void 
+Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount) const
+{
+  int nx = im.nx();
+  int ny = im.ny();
+  if (nx < 2 || ny < 2)
+      return;
+
+  int nsample = in_nsample;
+  if (nsample < 1)  
+    nsample = 1;
+
+  double dx = m_xmax - m_xmin;
+  double dy = m_ymax - m_ymin;
+  double xcent = m_xmin + dx / 2;
+  double ycent = m_ymin + dy / 2;
+  double phmlen = (dx > dy ? dx : dy);
+
+  double phmradius = phmlen / 2;
+
+  double xmin = xcent - phmradius;
+  double xmax = xcent + phmradius;
+  double ymin = ycent - phmradius;
+  double ymax = ycent + phmradius;
+
+  // Each pixel holds the average of the intensity of the cell with (ix,iy) at the center of the pixel
+  // Set major increments so that the last cell v[nx-1][ny-1] will start at xmax - xinc, ymax - yinc).
+  // Set minor increments so that sample points are centered in cell
+
+  double xinc = (xmax - xmin) / nx;
+  double yinc = (ymax - ymin) / ny;
+
+  double kxinc = xinc / nsample;               /* interval between samples */
+  double kyinc = yinc / nsample;
+  double kxofs = kxinc / 2;            /* offset of 1st point */
+  double kyofs = kyinc / 2;
+
+  im.setAxisExtent (xmin, xmax, ymin, ymax);
+  im.setAxisIncrement (xinc, yinc);
+
+  ImageFileArray v = im.getArray();
+
+  double x_start = xmin + (colStart * xinc);
+  for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
+    double x, y, xi, yi;
+    int ix, iy, kx, ky;
+    for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) {
+      for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
+       for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
+         for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
+           if ((*pelem)->isPointInside (xi, yi, PHM_COORD) == TRUE)
+             v[ix][iy] += (*pelem)->atten();
+       } // for kx
+      } /* for iy */
+    }  /* for ix */
+  }  /* for pelem */
+  
+
+  if (nsample > 1) {
+    double factor = 1.0 / (nsample * nsample);
+
+    for (int ix = 0; ix < colCount; ix++)
+      for (int iy = 0; iy < ny; iy++)
+       v[ix][iy] *= factor;
+  }
+}
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 // CLASS IDENTIFICATION
 //