r113: *** empty log message ***

[ctsim.git] / libctsim / phantom.cpp

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index 1508bea745912e3ffc07661c057f52c2096db4b1..e28074f5917f5555c0c017f9a62a4e6ef3c02d82 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.2 2000/06/19 19:54:23 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
@@ -315,6 +315,84 @@ 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 colStart, const int colCount, const int in_nsample, const int trace) const
+{
+  int nx = im.nx();
+  int ny = im.ny();
+  if (nx < 2 || ny < 2)
+      return;
+
+  if (nsample < 1)  
+        nsample = 1;
+
+  int nsample = in_nsample;
+  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
 //