Interim work on creating projections from an imagefile

[ctsim.git] / libctsim / phantom.cpp

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index ba4b0e43561de92191e50a6bf9f981cffdf41063..2b969e87e25128808cc399a2cd23bd95aa5d709a 100644 (file)
--- a/libctsim/phantom.cpp
+++ b/libctsim/phantom.cpp
@@ -78,6 +78,7 @@ Phantom::init ()
   m_composition = P_PELEMS;
   m_fail = false;
   m_id = PHM_INVALID;
+  m_im = NULL;
 }
 
 Phantom::~Phantom ()
@@ -85,9 +86,10 @@ Phantom::~Phantom ()
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
     delete *i;
   }
+  if (m_im)
+    delete m_im;
 }
 
-
 const char*
 Phantom::convertPhantomIDToName (int phmID)
 {
@@ -115,13 +117,13 @@ Phantom::convertNameToPhantomID (const char* const phmName)
 {
   int id = PHM_INVALID;
 
-  for (int i = 0; i < s_iPhantomCount; i++)
+  for (int i = 0; i < s_iPhantomCount; i++) {
     if (strcasecmp (phmName, s_aszPhantomName[i]) == 0) {
       id = i;
       break;
     }
-
-    return (id);
+  }
+  return (id);
 }
 
 
@@ -171,10 +173,10 @@ Phantom::createFromPhantom (const int phmid)
 
 
 /* METHOD IDENTIFICATION
-*   createFromFile          Add PhantomElements from file
+*   createFromPhmFile          Add PhantomElements from file
 *
 * SYNOPSIS
-*   createFromFile (filename)
+*   createFromPhmFile (filename)
 *
 * RETURNS
 *   true if pelem were added
@@ -182,7 +184,7 @@ Phantom::createFromPhantom (const int phmid)
 */
 
 bool
-Phantom::createFromFile (const char* const fname)
+Phantom::createFromPhmFile (const char* const fname)
 {
   bool bGoodFile = true;
   FILE *fp;
@@ -201,7 +203,7 @@ Phantom::createFromFile (const char* const fname)
     if (status == static_cast<int>(EOF))
       break;
     else if (status != 7) {
-      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromFile]", fname);
+      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromPhmFile]", fname);
       bGoodFile = false;
     }
     addPElem (pelemtype, cx, cy, u, v, rot, dens);
@@ -212,6 +214,32 @@ Phantom::createFromFile (const char* const fname)
   return (bGoodFile);
 }
 
+bool
+Phantom::createFromImageFile (const char* const fname)
+{
+  bool bGoodFile = true;
+
+  m_im = new ImageFile ();
+  if (! m_im || ! m_im->fileRead (fname)) {
+    sys_error (ERR_SEVERE, "Unable to read image file %s", fname);
+    bGoodFile = false;
+    delete m_im;
+    m_im = NULL;
+    m_fail = true;
+  } else {
+    m_name = fname;
+    m_id = -1;
+    m_fail = false;
+
+    if (! m_im->getAxisExtent(m_xmin, m_xmax, m_ymin, m_ymax)) {
+      m_xmax = m_im->nx() / 2.;  m_xmin = -m_xmax;
+      m_ymax = m_im->ny() / 2.;  m_ymin = -m_ymax;
+    }
+  }
+
+  return (bGoodFile);
+}
+
 bool
 Phantom::fileWrite (const char* const fname)
 {
@@ -507,6 +535,34 @@ Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const do
       *vCol++ = 0;
   }
 
+#if HAVE_OPENMP
+  double x_start = xmin + (colStart * xinc);
+  for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
+    const PhantomElement& rPElem = **pelem;
+    #pragma omp parallel for
+    for (int ix = 0; ix < colCount; ix++) {
+      double x = x_start + ix * xinc;
+      int iColStore = ix + iStorageOffset;
+      ImageFileColumn vCol = v[iColStore];
+
+      double y = ymin;
+      for (int iy = 0; iy < ny; iy++, y += yinc) {
+        double dAtten = 0;
+        double xi = x + kxofs;
+        for (int kx = 0; kx < nsample; kx++, xi += kxinc) {
+          double yi = y + kyofs;
+          for (int ky = 0; ky < nsample; ky++, yi += kyinc) {
+            if (rPElem.isPointInside (xi, yi, PHM_COORD))
+              dAtten += rPElem.atten();
+          } // ky
+        } // kx
+        *vCol++ += dAtten;
+      } /* iy */
+    }  /* ix */
+  }  /* pelem */
+
+#else
+
   double x_start = xmin + (colStart * xinc);
   for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
     const PhantomElement& rPElem = **pelem;
@@ -526,12 +582,13 @@ Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const do
       } /* for iy */
     }  /* for ix */
   }  /* for pelem */
-
+#endif
 
   if (nsample > 1) {
     double factor = 1.0 / static_cast<double>(nsample * nsample);
-
-
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (int ix = 0; ix < colCount; ix++) {
       int iColStore = ix + iStorageOffset;
       ImageFileColumn vCol = v[iColStore];