X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=e9e053d96cd0d19b357fa2efe2b5125ff465a07c;hp=7b13546bf19ab8c711e8ae1baa81f05fd8ab5b5c;hb=9776c9a12ba53419d34563a5ec57c90e3d6798f4;hpb=949b7c917af311cd38f911bc5d0af36f524a6086

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 7b13546..e9e053d 100644
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -6,9 +6,9 @@
 **   Date Started: Aug 84
 **
 **  This is part of the CTSim program
-**  Copyright (C) 1983-2000 Kevin Rosenberg
+**  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.6 2000/06/27 10:48:11 kevin Exp $
+**  $Id: projections.cpp,v 1.49 2001/02/22 18:22:40 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
@@ -26,33 +26,56 @@
 
 #include "ct.h"
 
+const kuint16 Projections::m_signature = ('P'*256 + 'J');
+
+const int Projections::POLAR_INTERP_INVALID = -1;
+const int Projections::POLAR_INTERP_NEAREST = 0;
+const int Projections::POLAR_INTERP_BILINEAR = 1;
+const int Projections::POLAR_INTERP_BICUBIC = 2;
+
+const char* const Projections::s_aszInterpName[] = 
+{
+  {"nearest"},
+  {"bilinear"},
+//  {"bicubic"},
+};
+
+const char* const Projections::s_aszInterpTitle[] = 
+{
+  {"Nearest"},
+  {"Bilinear"},
+//  {"Bicubic"},
+};
+
+const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*);
+
 
 /* NAME
- *    Projections		Constructor for projections matrix storage 
- *
- * SYNOPSIS
- *    proj = projections_create (filename, nView, nDet)
- *    Projections& proj		Allocated projections structure & matrix
- *    int nView			Number of rotated view
- *    int nDet			Number of detectors
- *
- */
+*    Projections		Constructor for projections matrix storage 
+*
+* SYNOPSIS
+*    proj = projections_create (filename, nView, nDet)
+*    Projections& proj		Allocated projections structure & matrix
+*    int nView			Number of rotated view
+*    int nDet			Number of detectors
+*
+*/
 
 Projections::Projections (const Scanner& scanner)
-  : m_projData(0)
+: m_projData(0)
 {
   initFromScanner (scanner);
 }
 
 
 Projections::Projections (const int nView, const int nDet)
-  : m_projData(0)
+: m_projData(0)
 {
   init (nView, nDet);
 }
 
 Projections::Projections (void)
-  : m_projData(0)
+: m_projData(0)
 {
   init (0, 0);
 }
@@ -62,27 +85,77 @@ Projections::~Projections (void)
   deleteProjData();
 }
 
+int
+Projections::convertInterpNameToID (const char* const interpName)
+{
+  int interpID = POLAR_INTERP_INVALID;
+  
+  for (int i = 0; i < s_iInterpCount; i++)
+    if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
+      interpID = i;
+      break;
+    }
+    
+    return (interpID);
+}
+
+const char*
+Projections::convertInterpIDToName (const int interpID)
+{
+  static const char *interpName = "";
+  
+  if (interpID >= 0 && interpID < s_iInterpCount)
+    return (s_aszInterpName[interpID]);
+  
+  return (interpName);
+}
+
+const char*
+Projections::convertInterpIDToTitle (const int interpID)
+{
+  static const char *interpTitle = "";
+  
+  if (interpID >= 0 && interpID < s_iInterpCount)
+    return (s_aszInterpTitle[interpID]);
+  
+  return (interpTitle);
+}
+
+
 
 void
 Projections::init (const int nView, const int nDet)
 {
+  m_label.setLabelType (Array2dFileLabel::L_HISTORY);
   m_nView = nView;
   m_nDet = nDet;
   newProjData ();
+  
+  time_t t = time (NULL);
+  tm* lt = localtime (&t);
+  m_year = lt->tm_year;
+  m_month = lt->tm_mon;
+  m_day = lt->tm_mday;
+  m_hour = lt->tm_hour;
+  m_minute = lt->tm_min;
+  m_second = lt->tm_sec;
 }
 
 void
 Projections::initFromScanner (const Scanner& scanner)
 {
+  m_label.setLabelType (Array2dFileLabel::L_HISTORY);
   deleteProjData();
   init (scanner.nView(), scanner.nDet());
-
-  m_phmLen = scanner.phmLen();
+  
   m_rotInc = scanner.rotInc();
   m_detInc = scanner.detInc();
+  m_geometry = scanner.geometry();
+  m_dFocalLength = scanner.focalLength();
+  m_dViewDiameter = scanner.viewDiameter();
   m_rotStart = 0;
-  m_detStart =  -scanner.radius() + (scanner.detInc() / 2);
-  m_phmLen = scanner.phmLen();
+  m_detStart =  -(scanner.detLen() / 2);
+  m_dFanBeamAngle = scanner.fanBeamAngle();
 }
 
 void
@@ -100,10 +173,10 @@ Projections::newProjData (void)
 {
   if (m_projData)
     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
-
+  
   if (m_nView > 0 && m_nDet) {
     m_projData = new DetectorArray* [m_nView];
-
+    
     for (int i = 0; i < m_nView; i++)
       m_projData[i] = new DetectorArray (m_nDet);
   }
@@ -111,12 +184,12 @@ Projections::newProjData (void)
 
 
 /* NAME
- *   projections_free			Free memory allocated to projections
- *
- * SYNOPSIS
- *   projections_free(proj)
- *   Projections& proj				Projectionss to be deallocated
- */
+*   projections_free			Free memory allocated to projections
+*
+* SYNOPSIS
+*   projections_free(proj)
+*   Projections& proj				Projectionss to be deallocated
+*/
 
 void 
 Projections::deleteProjData (void)
@@ -124,7 +197,7 @@ Projections::deleteProjData (void)
   if (m_projData != NULL) {
     for (int i = 0; i < m_nView; i++)
       delete m_projData[i];
-
+    
     delete m_projData;
     m_projData = NULL;
   }
@@ -132,9 +205,9 @@ Projections::deleteProjData (void)
 
 
 /* NAME
- *    Projections::headerWwrite         Write data header for projections file
- *
- */
+*    Projections::headerWwrite         Write data header for projections file
+*
+*/
 
 bool 
 Projections::headerWrite (fnetorderstream& fs)
@@ -151,18 +224,20 @@ Projections::headerWrite (fnetorderstream& fs)
   kuint16 _hour = m_hour;
   kuint16 _minute = m_minute;
   kuint16 _second = m_second;
-
+  
   kfloat64 _calcTime = m_calcTime;
   kfloat64 _rotStart = m_rotStart;
   kfloat64 _rotInc = m_rotInc;
   kfloat64 _detStart = m_detStart;
   kfloat64 _detInc = m_detInc;
-  kfloat64 _phmLen = m_phmLen;
-  
+  kfloat64 _viewDiameter = m_dViewDiameter;
+  kfloat64 _focalLength = m_dFocalLength;
+  kfloat64 _fanBeamAngle = m_dFanBeamAngle;
+
   fs.seekp(0);
   if (! fs)
     return false;
-
+  
   fs.writeInt16 (_hsize);
   fs.writeInt16 (_signature);
   fs.writeInt32 (_nView);
@@ -173,7 +248,9 @@ Projections::headerWrite (fnetorderstream& fs)
   fs.writeFloat64 (_rotInc);
   fs.writeFloat64 (_detStart);
   fs.writeFloat64 (_detInc);
-  fs.writeFloat64 (_phmLen);
+  fs.writeFloat64 (_viewDiameter);
+  fs.writeFloat64 (_focalLength);
+  fs.writeFloat64 (_fanBeamAngle);
   fs.writeInt16 (_year);
   fs.writeInt16 (_month);
   fs.writeInt16 (_day);
@@ -182,33 +259,32 @@ Projections::headerWrite (fnetorderstream& fs)
   fs.writeInt16 (_second);
   fs.writeInt16 (_remarksize);
   fs.write (m_remark.c_str(), _remarksize);
-
+  
   m_headerSize = fs.tellp();
   _hsize = m_headerSize;
   fs.seekp(0);
   fs.writeInt16 (_hsize);
   if (! fs)
-      return false;
+    return false;
   
   return true;
 }
 
 /* NAME
- *    projections_read_header         Read data header for projections file
- *
- */
+*    projections_read_header         Read data header for projections file
+*
+*/
 bool
 Projections::headerRead (fnetorderstream& fs)
 {
   kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
   kuint32 _nView, _nDet, _geom;
-  kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen;
+  kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _focalLength, _viewDiameter, _fanBeamAngle;
   
   fs.seekg(0);
   if (! fs)
-      return false;
-
-  off_t testPos;
+    return false;
+  
   fs.readInt16 (_hsize);
   fs.readInt16 (_signature);
   fs.readInt32 (_nView);
@@ -219,7 +295,9 @@ Projections::headerRead (fnetorderstream& fs)
   fs.readFloat64 (_rotInc);
   fs.readFloat64 (_detStart);
   fs.readFloat64 (_detInc);
-  fs.readFloat64 (_phmLen);
+  fs.readFloat64 (_viewDiameter);
+  fs.readFloat64 (_focalLength);
+  fs.readFloat64 (_fanBeamAngle);
   fs.readInt16 (_year);
   fs.readInt16 (_month);
   fs.readInt16 (_day);
@@ -227,30 +305,31 @@ Projections::headerRead (fnetorderstream& fs)
   fs.readInt16 (_minute);
   fs.readInt16 (_second);
   fs.readInt16 (_remarksize);
-
+  
   if (! fs) {
-      sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
-      return false;
+    sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
+    return false;
   }
-
+  
   if (_signature != m_signature) {
     sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
     return false;
   }
-
-  char remarkStorage[_remarksize+1];
-  fs.read (remarkStorage, _remarksize);
+  
+  char* pszRemarkStorage = new char [_remarksize+1];
+  fs.read (pszRemarkStorage, _remarksize);
   if (! fs) {
     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
     return false;
   }
-  remarkStorage[_remarksize] = 0;
-  m_remark = remarkStorage;
-
+  pszRemarkStorage[_remarksize] = 0;
+  m_remark = pszRemarkStorage;
+  delete pszRemarkStorage;
+  
   off_t _hsizeread = fs.tellg();
   if (!fs || _hsizeread != _hsize) {
-      sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize);
-      return false;
+    sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize);
+    return false;
   }
   
   m_headerSize = _hsize;
@@ -262,91 +341,200 @@ Projections::headerRead (fnetorderstream& fs)
   m_rotInc = _rotInc;
   m_detStart = _detStart;
   m_detInc = _detInc;
-  m_phmLen = _phmLen;
+  m_dFocalLength = _focalLength;
+  m_dViewDiameter = _viewDiameter;
+  m_dFanBeamAngle = _fanBeamAngle;
   m_year = _year;
   m_month = _month;
   m_day = _day;
   m_hour = _hour;
   m_minute = _minute;
   m_second = _second;
-
+  
   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
   m_label.setLabelString (m_remark);
   m_label.setCalcTime (m_calcTime);
   m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
-
+  
   return true;
 }
 
+bool
+Projections::read (const std::string& filename)
+{
+  return read (filename.c_str());
+}
+
+
 bool
 Projections::read (const char* filename)
 {
-  frnetorderstream fileRead (filename, ios::in | ios::binary);
   m_filename = filename;
-
-  if (! fileRead)
+#ifdef MSVC
+  frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary);
+#else
+  frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate);
+#endif
+  
+  if (fileRead.fail())
     return false;
-
+  
   if (! headerRead (fileRead))
     return false;
-
+  
   deleteProjData ();
   newProjData();
-
+  
   for (int i = 0; i < m_nView; i++) {
     if (! detarrayRead (fileRead, *m_projData[i], i))
       break;
   }
-
+  
   fileRead.close();
   return true;
 }
 
 
+bool 
+Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView)
+{
+  return copyViewData (filename.c_str(), os, startView, endView);
+}
+
+bool 
+Projections::copyViewData (const char* const filename, std::ostream& os, int startView, int endView)
+{
+  frnetorderstream is (filename, std::ios::in | std::ios::binary);
+  kuint16 sizeHeader, signature;
+  kuint32 _nView, _nDet;
+  
+  is.seekg (0);
+  if (is.fail()) {
+    sys_error (ERR_SEVERE, "Unable to read projection file %s", filename);
+    return false;
+  }
+
+  is.readInt16 (sizeHeader);
+  is.readInt16 (signature);
+  is.readInt32 (_nView);
+  is.readInt32 (_nDet);
+  int nView = _nView;
+  int nDet = _nDet;
+  
+  if (signature != m_signature) {
+    sys_error (ERR_SEVERE, "Illegal signature in projection file %s", filename);
+    return false;
+  }
+  
+  if (startView < 0)
+    startView = 0;
+  if (startView > nView - 1)
+    startView = nView;
+  if (endView < 0 || endView > nView - 1)
+    endView = nView - 1;
+  
+  if (startView > endView) { // swap if start > end
+    int tempView = endView;
+    endView = startView;
+    startView = tempView;
+  }
+  
+  int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
+  unsigned char* pViewData = new unsigned char [sizeView];
+  
+  for (int i = startView; i <= endView; i++) {
+    is.seekg (sizeHeader + i * sizeView);
+    is.read (reinterpret_cast<char*>(pViewData), sizeView);
+    os.write (reinterpret_cast<char*>(pViewData), sizeView);
+    if (is.fail() || os.fail())
+      break;
+  }
+  
+  delete pViewData;
+  if (is.fail()) 
+    sys_error (ERR_SEVERE, "Error reading projection file");
+  if (os.fail()) 
+    sys_error (ERR_SEVERE, "Error writing projection file");
+  
+  return (! (is.fail() | os.fail()));
+}
+
+bool 
+Projections::copyHeader (const std::string& filename, std::ostream& os)
+{
+  return copyHeader (filename.c_str(), os);
+}
+
+bool
+Projections::copyHeader (const char* const filename, std::ostream& os)
+{
+  frnetorderstream is (filename, std::ios::in | std::ios::binary);
+  kuint16 sizeHeader, signature;
+  is.readInt16 (sizeHeader);
+  is.readInt16 (signature);
+  is.seekg (0);
+  if (signature != m_signature) {
+    sys_error (ERR_SEVERE, "Illegal signature in projection file %s", filename);
+    return false;
+  }
+  
+  unsigned char* pHdrData = new unsigned char [sizeHeader];
+  is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
+  if (is.fail()) {
+    sys_error (ERR_SEVERE, "Error reading header");
+    return false;
+  }
+  
+  os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
+  if (os.fail()) {
+    sys_error (ERR_SEVERE, "Error writing header");
+    return false;
+  }
+  
+  return true;
+}
+
+bool
+Projections::write (const std::string& filename)
+{
+  return write (filename.c_str());
+}
+
 bool
 Projections::write (const char* filename)
 {
-  frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
+  frnetorderstream fs (filename, std::ios::out | std::ios::binary | std::ios::trunc | std::ios::ate);
   m_filename = filename;
   if (! fs) {
     sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
     return false;
   }
-
-  time_t t = time(NULL);
-  tm* lt = localtime(&t);
-  m_year = lt->tm_year;
-  m_month = lt->tm_mon;
-  m_day = lt->tm_mday;
-  m_hour = lt->tm_hour;
-  m_minute = lt->tm_min;
-  m_second = lt->tm_sec;
-
+  
   if (! headerWrite (fs))
-      return false;
-
+    return false;
+  
   if (m_projData != NULL) {
     for (int i = 0; i < m_nView; i++) {
       if (! detarrayWrite (fs, *m_projData[i], i))
-	break;
+        break;
     }
   }
   if (! fs)
     return false;
-
- fs.close();
-
+  
+  fs.close();
+  
   return true;
 }
 
 /* NAME
- *   detarrayRead		Read a Detector Array structure from the disk
- *
- * SYNOPSIS
- *   detarrayRead (proj, darray, view_num)
- *   DETARRAY *darray		Detector array storage location to be filled
- *   int      view_num		View number to read
- */
+*   detarrayRead		Read a Detector Array structure from the disk
+*
+* SYNOPSIS
+*   detarrayRead (proj, darray, view_num)
+*   DETARRAY *darray		Detector array storage location to be filled
+*   int      view_num		View number to read
+*/
 
 bool
 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
@@ -358,39 +546,39 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int
   DetectorValue* detval_ptr = darray.detValues();  
   kfloat64 view_angle;
   kuint32 nDet;
-
+  
   fs.seekg (start_data);
-
+  
   fs.readFloat64 (view_angle);
   fs.readInt32 (nDet);
   darray.setViewAngle (view_angle);
   //  darray.setNDet ( nDet);
   
-  for (int i = 0; i < nDet; i++) {
-      kfloat32 detval;
-      fs.readFloat32 (detval);
-      detval_ptr[i] = detval;
+  for (unsigned int i = 0; i < nDet; i++) {
+    kfloat32 detval;
+    fs.readFloat32 (detval);
+    detval_ptr[i] = detval;
   }
   if (! fs)
     return false;
-
+  
   return true;
 }
 
 
 /* NAME
- *   detarrayWrite			Write detector array data to the disk
- *
- * SYNOPSIS
- *   detarrayWrite (darray, view_num)
- *   DETARRAY *darray			Detector array data to be written
- *   int      view_num			View number to write
- *
- * DESCRIPTION
- *       This routine writes the detarray data from the disk sequentially to
- *    the file that was opened with open_projections().  Data is written in
- *    binary format.
- */
+*   detarrayWrite			Write detector array data to the disk
+*
+* SYNOPSIS
+*   detarrayWrite (darray, view_num)
+*   DETARRAY *darray			Detector array data to be written
+*   int      view_num			View number to write
+*
+* DESCRIPTION
+*       This routine writes the detarray data from the disk sequentially to
+*    the file that was opened with open_projections().  Data is written in
+*    binary format.
+*/
 
 bool
 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
@@ -408,205 +596,282 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co
     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
     return false;
   }
-
+  
   fs.writeFloat64 (view_angle);
   fs.writeInt32 (nDet);
-
-  for (int i = 0; i < nDet; i++) {
+  
+  for (unsigned int i = 0; i < nDet; i++) {
     kfloat32 detval = detval_ptr[i];
     fs.writeFloat32 (detval);
   }
-
+  
   if (! fs)
     return (false);
-
+  
   return true;
 }
 
 /* NAME
- *   prt_projections			Print projections data
- *
- * SYNOPSIS
- *   prt_projections (proj)
- *   Projections& proj			Projection data to be printed
- */
+*   printProjectionData			Print projections data
+*
+* SYNOPSIS
+*   printProjectionData ()
+*/
 
 void
-Projections::printProjectionData (void)
+Projections::printProjectionData ()
 {
-  printf("Projections Print\n\n");
-  printf("Description: %s\n", m_remark.c_str());
-  printf("nView = %d  nDet = %d\n", m_nView, m_nDet);
-  printf("rotStart = %8.4f   rotInc = %8.4f\n", m_rotStart, m_rotInc);
-  printf("detStart = %8.4f   detInc = %8.4f\n", m_detStart, m_detInc);
+  printProjectionData (0, nView() - 1);
+}
 
+void
+Projections::printProjectionData (int startView, int endView)
+{
+  printf("Projections Data\n\n");
+  printf("Description: %s\n", m_remark.c_str());
+  printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
+  printf("nView       = %8d           nDet = %8d\n", m_nView, m_nDet);
+  printf("focalLength = %8.4f ViewDiameter = %8.4f\n", m_dFocalLength, m_dViewDiameter);
+  printf("fanBeamAngle= %8.4f\n", convertRadiansToDegrees(m_dFanBeamAngle));
+  printf("rotStart    = %8.4f       rotInc = %8.4f\n", m_rotStart, m_rotInc);
+  printf("detStart    = %8.4f       detInc = %8.4f\n", m_detStart, m_detInc);
   if (m_projData != NULL) {
-      for (int ir = 0; ir < m_nView; ir++) {
-	  DetectorValue* detval = m_projData[ir]->detValues();
-	  for (int id = 0; id < m_projData[ir]->nDet(); id++)
-	      printf("%8.4f  ", detval[id]);
-	  printf("\n");
-      }
+    if (startView < 0)
+      startView = 0;
+    if (endView < 0)
+      endView = m_nView - 1;
+    if (startView > m_nView - 1)
+      startView = m_nView - 1;
+    if (endView > m_nView - 1)
+      endView = m_nView - 1;
+    for (int ir = startView; ir <= endView - 1; ir++) {
+      printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
+      DetectorValue* detval = m_projData[ir]->detValues();
+      for (int id = 0; id < m_projData[ir]->nDet(); id++)
+        printf("%8.4f  ", detval[id]);
+      printf("\n");
+    }
   }
 }
 
 void 
-Projections::printScanInfo (void) const
-{
-  printf ("Number of detectors: %d\n", m_nDet);
-  printf ("    Number of views: %d\n", m_nView);
-  printf ("             Remark: %s\n", m_remark.c_str());
-  printf ("             phmLen: %f\n", m_phmLen);
-  printf ("           detStart: %f\n", m_detStart);
-  printf ("             detInc: %f\n", m_detInc);
-  printf ("           rotStart: %f\n", m_rotStart);
-  printf ("             rotInc: %f\n", m_rotInc);
+Projections::printScanInfo (std::ostringstream& os) const
+{
+  os << "Number of detectors: " << m_nDet << "\n";
+  os << "Number of views: " << m_nView<< "\n";
+  os << "Description: " << m_remark.c_str()<< "\n";
+  os << "Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
+  os << "Focal Length: " << m_dFocalLength<< "\n";
+  os << "View Diameter: " << m_dViewDiameter<< "\n";
+  os << "Fan Beam Angle: " << convertRadiansToDegrees(m_dFanBeamAngle) << "\n";
+  os << "detStart: " << m_detStart<< "\n";
+  os << "detInc: " << m_detInc<< "\n";
+  os << "rotStart: " << m_rotStart<< "\n";
+  os << "rotInc: " << m_rotInc<< "\n";
 }
 
 
+bool 
+Projections::convertPolar (ImageFile& rIF, int iInterpolationID)
+{
+  unsigned int nx = rIF.nx();
+  unsigned int ny = rIF.ny();
+  ImageFileArray v = rIF.getArray();
+  ImageFileArray vImag = rIF.getImaginaryArray();
 
-/* NAME
- *   Projections::reconstruct      Reconstruct Image from Projections
- *
- * SYNOPSIS
- *   im = proj.reconstruct (im, filt_type, filt_param, interp_type)
- *   IMAGE *im			Output image
- *   int filt_type		Type of convolution filter to use
- *   double filt_param		Filter specific parameter
- *				Currently, used only with Hamming filters
- *   int interp_type		Type of interpolation method to use
- *
- * ALGORITHM
- *
- *	Calculate one-dimensional filter in spatial domain
- *	Allocate & clear (zero) the 2d output image array
- *      For each projection view
- *	    Convolve raysum array with filter
- *	    Backproject raysums and add (summate) to image array
- *      end
- */
-
-bool
-Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
-{
-  int nview = m_nView;
-  double det_inc = m_detInc;
-  double detlen = (m_nDet - 1) * det_inc;
-  int n_filtered_proj = m_nDet;
-  double filtered_proj [n_filtered_proj];   // convolved result
-
-#ifdef HAVE_BSPLINE_INTERP
-  int spline_order = 0, zoom_factor = 0;
-  if (interp_type == I_BSPLINE) {
-    zoom_factor = interp_param;
-    spline_order = 3;
-    zoom_factor = 3;
-    n_filtered_proj = (m_nDet - 1) * (zoom_factor + 1) + 1;
-  }
-#endif
+  if (! v || nx == 0 || ny == 0)
+    return false;
+  
+  Array2d<double> adView (nx, ny);
+  Array2d<double> adDet (nx, ny);
+  double** ppdView = adView.getArray();
+  double** ppdDet = adDet.getArray();
 
-  int n_vec_filter = 2 * m_nDet - 1;
-  double filterBW = 1. / det_inc;
-  double filterMin = -detlen;
-  double filterMax = detlen;
+  if (! calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet)) 
+    return false;
 
-  SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0);
-  if (filter.fail()) {
-      sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str());
-      return false;
+  std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+  unsigned int iView;
+  for (iView = 0; iView < m_nView; iView++) {
+    ppcDetValue[iView] = new std::complex<double> [m_nDet];
+    for (unsigned int iDet = 0; iDet < m_nDet; iDet++)
+      ppcDetValue[iView][iDet] = std::complex<double>(getDetectorArray (iView).detValues()[iDet], 0);
   }
 
-  if (trace)
-    cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
+  interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
 
-#if HAVE_SGP
-  SGP_ID gid;
-  double plot_xaxis [n_vec_filter]; 			// array for plotting 
+  for (iView = 0; iView < m_nView; iView++)
+    delete [] ppcDetValue[iView];
+  delete [] ppcDetValue;
 
-  if (trace > TRACE_TEXT)  {
-    int i;
-    double f;
-    double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
-    for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
-      plot_xaxis[i] = f;
-      
-    gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
-    cio_put_str ("Press any key to continue");
-    cio_kb_getc ();
-    sgp2_close (gid);
-  }
-  if (trace >= TRACE_TEXT) {
-    printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
-  }
-#endif  //HAVE_SGP
+  return true;
+}
 
-  Backprojector bj (*this, im, backprojectName, interpName);
-  if (bj.fail()) {
-    sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str());
+
+bool 
+Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
+{
+  unsigned int nx = rIF.nx();
+  unsigned int ny = rIF.ny();
+  ImageFileArray v = rIF.getArray();
+  if (! rIF.isComplex())
+    rIF.convertRealToComplex();
+  ImageFileArray vImag = rIF.getImaginaryArray();
+
+  if (! v || nx == 0 || ny == 0)
     return false;
+  
+#ifndef HAVE_FFT
+  return false;
+#else
+  Array2d<double> adView (nx, ny);
+  Array2d<double> adDet (nx, ny);
+  double** ppdView = adView.getArray();
+  double** ppdDet = adDet.getArray();
+
+  std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+  unsigned int iView;
+  double* pdDet = new double [m_nDet];
+  fftw_complex* pcIn = new fftw_complex [m_nDet];
+  fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
+
+  for (iView = 0; iView < m_nView; iView++) {
+    unsigned int iDet;
+    for (iDet = 0; iDet < m_nDet; iDet++) {
+      pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet];
+      pcIn[iDet].im = 0;
+    }
+    fftw_one (plan, pcIn, NULL);
+    ppcDetValue[iView] = new std::complex<double> [m_nDet];
+    for (iDet = 0; iDet < m_nDet; iDet++)
+      ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im); 
+    Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
   }
 
-  for (int iview = 0; iview < m_nView; iview++)  {
-    if (trace >= TRACE_TEXT) 
-      printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
-      
-    DetectorArray& darray = getDetectorArray (iview);
-    DetectorValue* detval = darray.detValues();
-
-    for (int j = 0; j < m_nDet; j++)
-      filtered_proj[j] = filter.convolve (detval, det_inc, j, m_nDet);
-
-#ifdef HAVE_SGP
-    if (trace >= TRACE_PLOT)  {
-      ezset  ("clear.");
-      ezset  ("xticks major 5.");
-      ezset  ("xlabel ");
-      ezset  ("ylabel ");
-      ezset  ("xlength .5.");
-      ezset  ("box.");
-      ezset  ("grid.");
-      ezset  ("ufinish yes.");
-      ezplot (detval, plot_xaxis, m_nDet);
-      ezset  ("clear.");
-      ezset  ("xticks major 5.");
-      ezset  ("xlabel ");
-      ezset  ("ylabel ");
-      ezset  ("ustart yes.");
-      ezset  ("xporigin .5.");
-      ezset  ("xlength .5.");
-      ezset ("box");
-      ezset ("grid");
-      gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
-    }
-#endif  //HAVE_SGP
+  fftw_destroy_plan (plan);  
+  delete [] pcIn;
+  
+  bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
 
-#ifdef HAVE_BSPLINE_INTERP
-    if (interp_type == I_BSPLINE) 
-	bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
-    
-#ifdef HAVE_SGP
-    if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
-	bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
-      ezplot_1d (filtered_proj, n_filtered_proj);
-    }
-#endif
+  if (! bError)
+    interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
+
+  for (iView = 0; iView < m_nView; iView++)
+    delete [] ppcDetValue[iView];
+  delete [] ppcDetValue;
+
+  return bError;
 #endif
+}
 
-    bj.BackprojectView (filtered_proj, darray.viewAngle());
 
-#ifdef HAVE_SGP
-    if (trace >= TRACE_PLOT) {
-      char str[256];
-      printf ("Do you want to exit with current pic (y/n) -- ");
-      fgets(str, sizeof(str), stdin);
-      sgp2_close (sgp2_get_active_win());
-      if (tolower(str[0]) == 'y') {
-	break;
-      }
-    } 
-#endif  //HAVE_SGP
+bool
+Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet)
+{
+  double xMin = -phmLen() / 2;
+  double xMax = xMin + phmLen();
+  double yMin = -phmLen() / 2;
+  double yMax = yMin + phmLen();
+  
+  double xInc = (xMax - xMin) / nx;	// size of cells
+  double yInc = (yMax - yMin) / ny;
+  
+  int iDetCenter = (m_nDet - 1) / 2;	// index refering to L=0 projection 
+
+  if (m_geometry != Scanner::GEOMETRY_PARALLEL) {
+    sys_error (ERR_WARNING, "convertPolar supports Parallel only");
+    return false;
+  }
+  
+  // Calculates polar coordinates (view#, det#) for each point on phantom grid
+  double x = xMin + xInc / 2;	// Rectang coords of center of pixel 
+  for (unsigned int ix = 0; ix < nx; x += xInc, ix++) {
+    double y = yMin + yInc / 2;
+    for (unsigned int iy = 0; iy < ny; y += yInc, iy++) {
+      double r = ::sqrt (x * x + y * y);
+      double phi = atan2 (y, x);
+
+      if (phi >= PI) {
+        phi -= PI;
+      } else if (phi < 0) {
+        phi += PI;
+      } else
+        r = -r;
+      
+      ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
+      ppdDet[ix][iy] = (r / m_detInc) + iDetCenter;
+    }
   }
 
   return true;
 }
 
+void
+Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
+     unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
+     double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID)
+{
+  for (unsigned int ix = 0; ix < ny; ix++) {
+    for (unsigned int iy = 0; iy < ny; iy++) {
+      if (iInterpolationID == POLAR_INTERP_NEAREST) {
+        int iView = nearest<int> (ppdView[ix][iy]);
+        int iDet = nearest<int> (ppdDet[ix][iy]);
+        if (iView == nView) {
+          iView = 0;
+       //   iDet = m_nDet - iDet;
+        }
+        if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) {
+          v[ix][iy] = ppcDetValue[iView][iDet].real();
+          if (vImag)
+            vImag[ix][iy] = ppcDetValue[iView][iDet].imag();
+        } else {
+          sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", 
+            ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+          v[ix][iy] = 0;
+        }
+      } else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
+        int iFloorView = static_cast<int>(ppdView[ix][iy]);
+        double dFracView = ppdView[ix][iy] - iFloorView;
+        int iFloorDet = static_cast<int>(ppdDet[ix][iy]);
+        double dFracDet = ppdDet[ix][iy] - iFloorDet;
+
+        if (iFloorDet >= 0 && iFloorView >= 0) { 
+          std::complex<double> v1 = ppcDetValue[iFloorView][iFloorDet];
+          std::complex<double> v2, v3, v4;
+          if (iFloorView < nView - 1)
+            v2 = ppcDetValue[iFloorView + 1][iFloorDet];
+          else 
+            v2 = ppcDetValue[0][iFloorDet];
+          if (iFloorDet < nDet - 1) 
+            v4 = ppcDetValue[iFloorView][iFloorDet+1];
+          else
+            v4 = v1;
+          if (iFloorView < nView - 1 && iFloorDet < nDet - 1)
+            v3 = ppcDetValue [iFloorView+1][iFloorDet+1];
+          else if (iFloorView < nView - 1)
+            v3 = v2;
+          else 
+            v3 = ppcDetValue[0][iFloorDet+1];
+          std::complex<double> vInterp = (1 - dFracView) * (1 - dFracDet) * v1 +
+            dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 +
+            dFracDet * (1 - dFracView) * v4;
+          v[ix][iy] = vInterp.real();
+          if (vImag)
+            vImag[ix][iy] = vInterp.imag();
+        } else {
+          sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", 
+            ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+          v[ix][iy] = 0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+        }
+      } else if (iInterpolationID == POLAR_INTERP_BICUBIC) {
+        v[ix][iy] =0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+      }
+    }
+  }
+}
+
+