changed sinogram from std::vector to Eigen matrix

include/Sinogram.h

@@ -38,7 +38,7 @@ class Sinogram
     Eigen::Matrix<uint32_t, 28, 1> dilationOffsets; // actually 27, but padded to be fixed-size-vectorizeable
     Eigen::Array<float, 1, Eigen::Dynamic> sinah, cosah;
     float realToMatrixScaling, realToMatrixOffset;
-    std::vector<uint8_t> sinogram, sinogram_oneMeasuredPeak;
+    Eigen::Matrix<uint8_t, Eigen::Dynamic, 1> sinogram, sinogram_oneMeasuredPeak;
 
   public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW

src/Sinogram.cpp

@@ -39,8 +39,8 @@ void Sinogram::setSinogramAngleResolution(float angleResolution_deg)
         throw BadInputException("sinogram too big to fit in uint32 adress range!");
     }
 
-    sinogram.resize(sinogramElementsCount);
-    sinogram_oneMeasuredPeak.resize(sinogramElementsCount);
+    sinogram.resize(sinogramElementsCount, 1);
+    sinogram_oneMeasuredPeak.resize(sinogramElementsCount, 1);
 
     sinogramCenter = sinogramSize / 2;
     sinogramScale = 1 / ((float)sinogramSize_exact / 2);
@@ -70,14 +70,11 @@ void Sinogram::setSinogramAngleResolution(float angleResolution_deg)
 
 void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms)
 {
-    fill(sinogram.begin(), sinogram.end(), 0);
+    sinogram.setZero();
 
     EigenSTL::vector_Matrix3Xf candidateReflectionsDirections;
     reflectionsInRangeFinder.getReflectionsInRanges(candidateReflectionsDirections, ucsBorderNorms);
 
-    Map<Matrix<uint8_t, Dynamic, 1>> sinogram_matrixMapped(sinogram.data(), sinogram.size()),
-        sinogram_oneMeasuredPeak_matrixMapped(sinogram_oneMeasuredPeak.data(), sinogram_oneMeasuredPeak.size());
-
     Matrix3Xf n(3, anglesCount);
     Array<float, 1, Eigen::Dynamic> gammah(anglesCount); // gamma half
     Matrix<uint32_t, 3, Eigen::Dynamic> validSinogramPoints_matrix(3, anglesCount);
@@ -100,7 +97,7 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige
         if (candidateReflectionDirections.cols() == 0)
             continue;
 
-        fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
+        sinogram_oneMeasuredPeak.setZero();
         int candidatePeaksCount = candidateReflectionDirections.cols();
         for (int candidatePeakNumber = 0; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber++)
         {
@@ -129,7 +126,7 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige
             }
         }
 
-        sinogram_matrixMapped += sinogram_oneMeasuredPeak_matrixMapped;
+        sinogram += sinogram_oneMeasuredPeak;
     }
 }
 
@@ -186,14 +183,11 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, co
     vector<thread> threads;
     threads.reserve(slaveThreadCount);
 
-    fill(sinogram.begin(), sinogram.end(), 0);
+    sinogram.setZero();
 
     EigenSTL::vector_Matrix3Xf candidateReflectionsDirections;
     reflectionsInRangeFinder.getReflectionsInRanges(candidateReflectionsDirections, ucsBorderNorms);
 
-    Map<Matrix<uint8_t, Dynamic, 1>> sinogram_matrixMapped(sinogram.data(), sinogram.size()),
-        sinogram_oneMeasuredPeak_matrixMapped(sinogram_oneMeasuredPeak.data(), sinogram_oneMeasuredPeak.size());
-
     int measuredPeaksCount = ucsDirections.cols();
     cout << "peak count used for indexing: " << measuredPeaksCount << endl;
     for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
@@ -209,7 +203,7 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, co
         if (candidateReflectionDirections.cols() == 0)
             continue;
 
-        fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
+        sinogram_oneMeasuredPeak.setZero();
 
         for (int threadNumber = 0; threadNumber < slaveThreadCount; threadNumber++)
         {
@@ -221,7 +215,7 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, co
         }
         threads.clear();
 
-        sinogram_matrixMapped += sinogram_oneMeasuredPeak_matrixMapped;
+        sinogram += sinogram_oneMeasuredPeak;
     }
 }
 
@@ -232,15 +226,14 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c
     vector<thread> threads;
     threads.reserve(slaveThreadCount);
 
-    vector<uint8_t> sinogram_oneMeasuredPeak_local(sinogram_oneMeasuredPeak.size(), 0);
+    Matrix<uint8_t, Eigen::Dynamic, 1> sinogram_oneMeasuredPeak_local(sinogram_oneMeasuredPeak.size(), 1);
+    sinogram_oneMeasuredPeak_local.setZero();
 
-    fill(sinogram.begin(), sinogram.end(), 0);
+    sinogram.setZero();
 
     EigenSTL::vector_Matrix3Xf candidateReflectionsDirections;
     reflectionsInRangeFinder.getReflectionsInRanges(candidateReflectionsDirections, ucsBorderNorms);
 
-    Map<Matrix<uint8_t, Dynamic, 1>> sinogram_matrixMapped(sinogram.data(), sinogram.size());
-
     int measuredPeaksCount = ucsDirections.cols();
     cout << "peak count used for indexing: " << measuredPeaksCount << endl;
     for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
@@ -250,15 +243,13 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c
             cout << "computed " << 100 * (float)measuredPeakNumber / measuredPeaksCount << "%" << endl;
         }
 
-
         Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis
 
         Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber];
         if (candidateReflectionDirections.cols() == 0)
             continue;
 
-        fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
-
+        sinogram_oneMeasuredPeak.setZero();
 
         for (int threadNumber = 0; threadNumber < slaveThreadCount; threadNumber++)
         {
@@ -267,11 +258,8 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c
 
         if (measuredPeakNumber > 0)
         {
-            Map<Matrix<uint8_t, Dynamic, 1>> sinogram_oneMeasuredPeakLocal_matrixMapped(sinogram_oneMeasuredPeak_local.data(),
-                                                                                        sinogram_oneMeasuredPeak_local.size());
-
-            sinogram_matrixMapped += sinogram_oneMeasuredPeakLocal_matrixMapped;
-            fill(sinogram_oneMeasuredPeak_local.begin(), sinogram_oneMeasuredPeak_local.end(), 0);
+            sinogram += sinogram_oneMeasuredPeak_local;
+            sinogram_oneMeasuredPeak_local.setZero();
         }
 
         for (int threadNumber = 0; threadNumber < slaveThreadCount; threadNumber++)
@@ -282,13 +270,14 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c
 
         sinogram_oneMeasuredPeak_local.swap(sinogram_oneMeasuredPeak);
     }
-    Map<Matrix<uint8_t, Dynamic, 1>> sinogram_oneMeasuredPeakLocal_matrixMapped(sinogram_oneMeasuredPeak_local.data(), sinogram_oneMeasuredPeak_local.size());
-    sinogram_matrixMapped += sinogram_oneMeasuredPeakLocal_matrixMapped;
+    sinogram += sinogram_oneMeasuredPeak_local;
 }
 
 void Sinogram::getBestRotation(AngleAxisf& bestRotation)
 {
-    uint32_t maxElementIdx = distance(sinogram.begin(), max_element(sinogram.begin(), sinogram.end()));
+    uint32_t maxElementIdx, dummy;
+    sinogram.maxCoeff(&maxElementIdx, &dummy);
+
     Matrix<uint32_t, 3, 1> maxElementSub;
     maxElementSub(2) = maxElementIdx / strides(1);
     maxElementIdx -= maxElementSub(2) * strides(1);