Interim work on creating projections from an imagefile

[ctsim.git] / libctsim / phantom.cpp

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index f02877775729f2f1a4781e78eec84ad2d9f1ef56..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)
 {
@@ -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)
 {
@@ -516,13 +544,14 @@ Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const do
       double x = x_start + ix * xinc;
       int iColStore = ix + iStorageOffset;
       ImageFileColumn vCol = v[iColStore];
-      for (int iy = 0; iy < ny; iy++) {
-        double y = ymin + iy * yinc;
+
+      double y = ymin;
+      for (int iy = 0; iy < ny; iy++, y += yinc) {
         double dAtten = 0;
-        for (int kx = 0; kx < nsample; kx++) {
-          double xi = x + kxofs + kxinc * kx;
-          for (int ky = 0; ky < nsample; ky++) {
-            double yi = y + kyofs + ky * kyinc;
+        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
@@ -531,9 +560,9 @@ Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const do
       } /* 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;
@@ -559,7 +588,7 @@ Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const do
     double factor = 1.0 / static_cast<double>(nsample * nsample);
 #if HAVE_OPENMP
     #pragma omp parallel for
-#endif    
+#endif
     for (int ix = 0; ix < colCount; ix++) {
       int iColStore = ix + iStorageOffset;
       ImageFileColumn vCol = v[iColStore];