projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r7061: initial property settings
[ctsim.git] / libctsim / projections.cpp
diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index b2d3c4aabc39cfcc37cc1eea20887322a5c9eebf..871f62e6e52a4543df9262355872d3b22957f392 100644 (file)
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.76 2001/09/24 11:29:41 kevin Exp $
+**  $Id$
 **
 **  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
@@ -25,6 +25,8 @@
 ******************************************************************************/
 
 #include "ct.h"
+#include <ctime>
+#include "interpolator.h"
 
 const kuint16 Projections::m_signature = ('P'*256 + 'J');
 
@@ -35,15 +37,15 @@ const int Projections::POLAR_INTERP_BICUBIC = 2;
 
 const char* const Projections::s_aszInterpName[] = 
 {
-  {"nearest"},
-  {"bilinear"},
+  "nearest",
+  "bilinear",
 //  {"bicubic"},
 };
 
 const char* const Projections::s_aszInterpTitle[] = 
 {
-  {"Nearest"},
-  {"Bilinear"},
+  "Nearest",
+  "Bilinear",
 //  {"Bicubic"},
 };
 
@@ -650,7 +652,7 @@ Projections::read (const char* filename)
 #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);
+  frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary); // | std::ios::nocreate);
 #endif
   
   if (fileRead.fail())
@@ -970,11 +972,11 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID)
   double** ppdDet = adDet.getArray();
 
   std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
-  unsigned int iView;
+  int iView;
   for (iView = 0; iView < pProj->m_nView; iView++) {
     ppcDetValue[iView] = new std::complex<double> [pProj->m_nDet];
     DetectorValue* detval = pProj->getDetectorArray (iView).detValues();
-    for (unsigned int iDet = 0; iDet < pProj->m_nDet; iDet++)
+    for (int iDet = 0; iDet < pProj->m_nDet; iDet++)
       ppcDetValue[iView][iDet] = std::complex<double>(detval[iDet], 0);
   }
 
@@ -1015,45 +1017,45 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
   if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR)
     pProj = interpolateToParallel();
 
-  int iInterpDet = nx;
-//  int iInterpDet = pProj->m_nDet;
+  int iInterpDet = static_cast<int>(static_cast<double>(sqrt(nx*nx+ny*ny)));
   int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
-
+  double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.5);
   double dZeropadRatio = static_cast<double>(iNumInterpDetWithZeros) / static_cast<double>(iInterpDet);
 
   fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
 
   fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros];
   std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
-  double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
+  //double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1);
+  double dInterpScale = pProj->m_nDet / static_cast<double>(iInterpDet);
   
   double dFFTScale = 1. / static_cast<double>(iInterpDet * iInterpDet);
   int iMidPoint = iInterpDet / 2;
   double dMidPoint = static_cast<double>(iInterpDet) / 2.;
   int iZerosAdded = iNumInterpDetWithZeros - iInterpDet;
 
-  // For each view, interpolate to nx length, shift to center at origin, and FFt transform
-  for (unsigned int iView = 0; iView < m_nView; iView++) {
+  // For each view, interpolate, shift to center at origin, and FFT
+  for (int iView = 0; iView < m_nView; iView++) {
     DetectorValue* detval = pProj->getDetectorArray(iView).detValues();
     LinearInterpolator<DetectorValue> projInterp (detval, pProj->m_nDet);
-    for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) {
+    for (int iDet = 0; iDet < iInterpDet; iDet++) {
       double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale;
-      pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale;
+      pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dProjScale;
       pcIn[iDet].im = 0;
     }
 
     Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet);
     if (iZerosAdded > 0) {
-      for (unsigned int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++)
+       for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--)
         pcIn[iDet1+iZerosAdded] = pcIn[iDet1];
-      for (unsigned int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) 
+      for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) 
         pcIn[iDet2].re = pcIn[iDet2].im = 0;
     }
 
     fftw_one (plan, pcIn, NULL);
 
     ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
-    for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
+    for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
       ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); 
     }
 
@@ -1090,7 +1092,6 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
                                         int iNumDetWithZeros, double dZeropadRatio, double dDetInc)
 {
   double dLength = viewDiameter();
-//  double dLength = phmLen();
   double xMin = -dLength / 2;
   double xMax = xMin + dLength;
   double yMin = -dLength / 2;
@@ -1106,10 +1107,10 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
   double xInc = (xMax - xMin) / nx;    // size of cells
   double yInc = (yMax - yMin) / ny;
 
+  double dDetCenter = (iNumDetWithZeros - 1) / 2.;     // index refering to L=0 projection 
   // +1 is correct for frequency data, ndet-1 is correct for projections
-  int iDetCenter = (iNumDetWithZeros - 1) / 2; // index refering to L=0 projection 
-  if (isEven (iNumDetWithZeros))
-    iDetCenter = (iNumDetWithZeros + 1) / 2;   
+  //  if (isEven (iNumDetWithZeros))
+  //    dDetCenter = (iNumDetWithZeros + 0) / 2;       
 
   // Calculates polar coordinates (view#, det#) for each point on phantom grid
   double x = xMin + xInc / 2;  // Rectang coords of center of pixel 
@@ -1119,15 +1120,15 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
       double r = ::sqrt (x * x + y * y);
       double phi = atan2 (y, x);
 
-      if (phi < 0)
+      if (phi <= -m_rotInc / 2)
         phi += TWOPI;
-      if (phi >= PI) {
+      if (phi >= PI - (m_rotInc / 2)) {
         phi -= PI;
         r = -r;
       }
       
       ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
-      ppdDet[ix][iy] = (r / dDetInc) + iDetCenter;
+      ppdDet[ix][iy] = (r / dDetInc) + dDetCenter;
     }
   }
 }
@@ -1139,30 +1140,29 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
 {
   typedef std::complex<double> complexValue;
 
-  BilinearInterpolator<complexValue>* pBilinear;  
+  BilinearPolarInterpolator<complexValue>* pBilinear = NULL;  
+  BicubicPolyInterpolator<complexValue>* pBicubic = NULL;  
   if (iInterpolationID == POLAR_INTERP_BILINEAR)
-    pBilinear = new BilinearInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
-
-  BicubicPolyInterpolator<complexValue>* pBicubic;  
-  if (iInterpolationID == POLAR_INTERP_BICUBIC)
+    pBilinear = new BilinearPolarInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
+  else if (iInterpolationID == POLAR_INTERP_BICUBIC)
     pBicubic = new BicubicPolyInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
-
+  
   for (unsigned int ix = 0; ix < ny; ix++) {
     for (unsigned int iy = 0; iy < ny; iy++) {
-
       if (iInterpolationID == POLAR_INTERP_NEAREST) {
         unsigned int iView = nearest<int> (ppdView[ix][iy]);
         unsigned int iDet = nearest<int> (ppdDet[ix][iy]);
-        if (iView == nView) {
-          iView = 0;
-          iDet = m_nDet - iDet;
-        }
+        if (iView == nView)
+         iView = 0;
         if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) {
           v[ix][iy] = ppcDetValue[iView][iDet].real();
           if (vImag)
             vImag[ix][iy] = ppcDetValue[iView][iDet].imag();
-        } else
+        } else {
           v[ix][iy] = 0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+       }
 
       } else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
         std::complex<double> vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]);
@@ -1206,22 +1206,22 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa
   for (int iv = 0; iv < iNViews; iv++) {
     unsigned char* pArgBase = pData + lDataPos;
     unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p);
-    long lProjNumber = *reinterpret_cast<long*>(p);
+    // long lProjNumber = *reinterpret_cast<long*>(p);
 
     p = pArgBase+20;  SwapBytes4IfLittleEndian (p);
     long lEscale = *reinterpret_cast<long*>(p);
 
     p = pArgBase+28;  SwapBytes4IfLittleEndian (p);
-    long lTime = *reinterpret_cast<long*>(p);
+    // long lTime = *reinterpret_cast<long*>(p);
 
     p = pArgBase + 4; SwapBytes4IfLittleEndian (p);
     double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
 
     p = pArgBase+12; SwapBytes4IfLittleEndian (p);
-    double dAlign = *reinterpret_cast<float*>(p);
+    // double dAlign = *reinterpret_cast<float*>(p);
 
     p = pArgBase + 16; SwapBytes4IfLittleEndian (p);
-    double dMaxValue = *reinterpret_cast<float*>(p);
+    // double dMaxValue = *reinterpret_cast<float*>(p);
 
     DetectorArray& detArray = getDetectorArray (iv);
     detArray.setViewAngle (dAlpha);
@@ -1292,9 +1292,9 @@ Projections::interpolateToParallel () const
 
     double dViewAngle = m_rotStart;
     int iLastFloor = -1;
+    LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
     for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
       DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
-      LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
       detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor);
     }
   }
@@ -1337,8 +1337,8 @@ Projections::interpolateToParallel () const
 ///////////////////////////////////////////////////////////////////////////////
 
 ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange)
-: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
-  m_iThetaRange (iThetaRange), m_pCoordinates(NULL)
+: m_pCoordinates(NULL), m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
+  m_iThetaRange (iThetaRange)
 {
   int iGeometry = pProjections->geometry();
   double dDetInc = pProjections->detInc();
Computed Tomography simulator and reconstructor