X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=c55aa37a5fcc94b6292ac9b912f93e858bbb5b7f;hp=d6a484335943945d8f4ffa02a6bb77e495412fb7;hb=999a754d1519a49ca062ee87b22bf601c1ee9f21;hpb=d850a3477e9ccaecfa85e00bc619848fcc29bdb6

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index d6a4843..c55aa37 100644
--- 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.56 2001/03/10 23:56:58 kevin Exp $
+**  $Id: projections.cpp,v 1.57 2001/03/11 06:34:37 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
@@ -150,12 +150,12 @@ Projections::initFromScanner (const Scanner& scanner)
   
   m_rotInc = scanner.rotInc();
   m_detInc = scanner.detInc();
+  m_detStart =  scanner.detStart();
   m_geometry = scanner.geometry();
   m_dFocalLength = scanner.focalLength();
   m_dSourceDetectorLength = scanner.sourceDetectorLength();
   m_dViewDiameter = scanner.viewDiameter();
   m_rotStart = 0;
-  m_detStart =  -(scanner.detLen() / 2);
   m_dFanBeamAngle = scanner.fanBeamAngle();
 }
 
@@ -968,34 +968,114 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa
   return true;
 }
 
+Projections*
+Projections::interpolateToParallel ()
+{
+  if (m_geometry == Scanner::GEOMETRY_PARALLEL)
+    return this;
+
+  int nDet = m_nDet;
+  int nView = m_nView;
+  Projections* pProjNew = new Projections (nView, nDet);
+  pProjNew->m_geometry = Scanner::GEOMETRY_PARALLEL;
+  pProjNew->m_dFocalLength = m_dFocalLength;
+  pProjNew->m_dSourceDetectorLength = m_dSourceDetectorLength;
+  pProjNew->m_dViewDiameter = m_dViewDiameter;
+  pProjNew->m_dFanBeamAngle = m_dFanBeamAngle;
+  pProjNew->m_calcTime  = 0;
+  pProjNew->m_remark = m_remark;
+  pProjNew->m_remark += "; Interpolate to Parallel";
+  pProjNew->m_label.setLabelType (Array2dFileLabel::L_HISTORY);
+  pProjNew->m_label.setLabelString (pProjNew->m_remark);
+  pProjNew->m_label.setCalcTime (pProjNew->m_calcTime);
+  pProjNew->m_label.setDateTime (pProjNew->m_year, pProjNew->m_month, pProjNew->m_day, pProjNew->m_hour, pProjNew->m_minute, pProjNew->m_second);
+
+  pProjNew->m_rotStart = 0;
+  pProjNew->m_rotInc = PI / nView;;
+  pProjNew->m_detStart = -m_dViewDiameter / 2;
+  pProjNew->m_detInc = m_dViewDiameter / nDet;
+  if (nDet % 2 == 0) // even
+    pProjNew->m_detInc = m_dViewDiameter / (nDet - 1);
+
+  ParallelRaysums parallel (this);
+
+  double* pdThetaValuesForT = new double [pProjNew->nView()];
+  double* pdRaysumsForT = new double [pProjNew->nView()];
+
+  // interpolate to evenly spaced theta (views)
+  double dDetPos = pProjNew->m_detStart;
+  for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
+    parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT);
+
+    double dViewAngle = m_rotStart;
+    for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
+      DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
+
+      detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle);
+    }
+  }
+  delete pdThetaValuesForT;
+  delete pdRaysumsForT;
+
+  // interpolate to evenly space t (detectors)
+  double* pdOriginalDetPositions = new double [pProjNew->nDet()];
+  parallel.getDetPositions (pdOriginalDetPositions);
+
+  double* pdDetValueCopy = new double [pProjNew->nDet()];
+  double dViewAngle = m_rotStart;
+  for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
+    DetectorArray& detArray = pProjNew->getDetectorArray (iV);
+    DetectorValue* detValues = detArray.detValues();
+    detArray.setViewAngle (dViewAngle);
+
+    for (int i = 0; i < pProjNew->nDet(); i++)
+      pdDetValueCopy[i] = detValues[i];
+
+    double dDetPos = pProjNew->m_detStart;
+    for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
+      detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos);
+    }
+  }
+  delete pdDetValueCopy;
+  delete pdOriginalDetPositions;
+
+  return pProjNew;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Class ParallelRaysums
+//
+// Used for converting divergent beam raysums into Parallel raysums
+//
+///////////////////////////////////////////////////////////////////////////////
 
 ParallelRaysums::ParallelRaysums (Projections* pProjections)
-: m_iNumCoordinates(0)
+: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet())
 {
-  int nDet = pProjections->nDet();
-  int nView = pProjections->nView();
   int iGeometry = pProjections->geometry();
   double dDetInc = pProjections->detInc();
   double dDetStart = pProjections->detStart();
   double dFocalLength = pProjections->focalLength();
 
-  m_iNumCoordinates =  nDet * nView;
+  m_iNumCoordinates =  m_iNumView * m_iNumDet;
   m_vecpCoordinates.reserve (m_iNumCoordinates);
   for (int i = 0; i < m_iNumCoordinates; i++)
     m_vecpCoordinates[i] = new ParallelRaysumCoordinate;
 
   int iCoordinate = 0;
-  for (int iV = 0; iV < nView; iV++) {
+  for (int iV = 0; iV < m_iNumView; iV++) {
     double dViewAngle = pProjections->getDetectorArray(iV).viewAngle();
+    const DetectorValue* detValues = pProjections->getDetectorArray(iV).detValues();
 
     double dDetPos = dDetStart;
-    for (int iD = 0; iD < nDet; iD++) {
+    for (int iD = 0; iD < m_iNumDet; iD++) {
       ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
 
       if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
         pC->m_dTheta = normalizeAngle (dViewAngle);
         pC->m_dT = dDetPos;
-
       } else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
         double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
         pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle);
@@ -1006,6 +1086,11 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections)
         pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos);
         pC->m_dT = dFocalLength * sin (dDetPos);        
       }
+      if (pC->m_dTheta >= PI) {  // convert T/Theta to 0-PI interval
+        pC->m_dTheta -= PI;
+        pC->m_dT = -pC->m_dT - pProjections->detInc();
+      }
+      pC->m_dRaysum = detValues[iD];
       dDetPos += dDetInc;
     }
   }
@@ -1021,7 +1106,7 @@ ParallelRaysums::CoordinateContainer&
 ParallelRaysums::getSortedByTheta()
 {
   if (m_vecpSortedByTheta.size() == 0) {
-    m_vecpSortedByTheta.reserve (m_iNumCoordinates);
+    m_vecpSortedByTheta.resize (m_iNumCoordinates);
     for (int i = 0; i < m_iNumCoordinates; i++)
       m_vecpSortedByTheta[i] = m_vecpCoordinates[i];
     std::sort (m_vecpSortedByTheta.begin(), m_vecpSortedByTheta.end(), ParallelRaysumCoordinate::compareByTheta);
@@ -1034,7 +1119,7 @@ ParallelRaysums::CoordinateContainer&
 ParallelRaysums::getSortedByT()
 {
   if (m_vecpSortedByT.size() == 0) {
-    m_vecpSortedByT.reserve (m_iNumCoordinates);
+    m_vecpSortedByT.resize (m_iNumCoordinates);
     for (int i = 0; i < m_iNumCoordinates; i++)
       m_vecpSortedByT[i] = m_vecpCoordinates[i];
     std::sort (m_vecpSortedByT.begin(), m_vecpSortedByT.end(), ParallelRaysumCoordinate::compareByT);
@@ -1068,3 +1153,56 @@ ParallelRaysums::getLimits (double* dMinT, double* dMaxT, double* dMinTheta, dou
       *dMaxTheta = dTheta;
   }
 }
+
+void
+ParallelRaysums::getThetaAndRaysumsForT (int iTheta, double* pTheta, double* pRaysum)
+{
+  const CoordinateContainer& coordsT = getSortedByT();
+
+  int iBase = iTheta * m_iNumView;
+  for (int i = 0; i < m_iNumView; i++) {
+    int iPos = iBase + i;
+    pTheta[i] = coordsT[iPos]->m_dTheta;
+    pRaysum[i] = coordsT[iPos]->m_dRaysum;
+  }
+}
+
+void
+ParallelRaysums::getDetPositions (double* pdDetPos)
+{
+  const CoordinateContainer& coordsT = getSortedByT();
+
+  int iPos = 0;
+  for (int i = 0; i < m_iNumDet; i++) {
+    pdDetPos[i] = coordsT[iPos]->m_dT;
+    iPos += m_iNumView;
+  }
+}
+
+// locate by bisection, O(log2(n))
+double
+ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX)
+{
+  int iLower = -1;
+  int iUpper = n;
+
+  while (iUpper - iLower > 1) {
+    int iMiddle = (iUpper + iLower) >> 1;
+    if (dX >= pdX[iMiddle])
+      iLower = iMiddle;
+    else
+      iUpper = iMiddle;
+  }
+  if (dX <= pdX[0])
+    return pdY[0];
+  else if (dX >= pdX[n-1])
+    return pdY[1];
+
+  if (iLower < 0 || iLower >= n) {
+    sys_error (ERR_SEVERE, "Coordinate out of range [locateThetaBase]");
+    return 0;
+  }
+
+  return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower]));
+}
+