diff --git a/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVectorTest.cpp b/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVectorTest.cpp
index 4eb6336ccf..90fcc05897 100644
--- a/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVectorTest.cpp
+++ b/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVectorTest.cpp
@@ -1,386 +1,251 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkTestFixture.h>
 #include <mitkTestingMacros.h>
+#include <math.h>
 
 #include "mitkPhotoacousticSmartVector.h"
 
 class mitkPhotoacousticVectorTestSuite : public mitk::TestFixture
 {
   CPPUNIT_TEST_SUITE(mitkPhotoacousticVectorTestSuite);
 
   MITK_TEST(TestNormalizedVector);
   MITK_TEST(TestRotatedVectorZeroDegrees);
   MITK_TEST(TestRotatedVectorPositiveDegrees);
   MITK_TEST(TestRotatedVectorZeroDegrees);
   MITK_TEST(TestRotatedVectorAroundAllAxis);
   MITK_TEST(TestScaledVector);
   MITK_TEST(TestClonedVector);
 
   CPPUNIT_TEST_SUITE_END();
 
 private:
 
   mitk::PhotoacousticSmartVector::Pointer m_TestVector;
   mitk::PhotoacousticSmartVector::Pointer m_TestReturnVector;
 
-  double DIF_VAL;
+  const double DIF_VAL = 0.001;
+  const double TWO_PI = 6.283185;
 
 public:
 
   void setUp()
   {
     m_TestVector = mitk::PhotoacousticSmartVector::New();
     m_TestReturnVector = mitk::PhotoacousticSmartVector::New();
   }
 
   void TestNormalizedVector()
   {
     std::stringstream output;
     int a = 2;
     int b = 3;
     int c = 4;
 
     m_TestVector->SetElement(0,a);
     m_TestVector->SetElement(1,b);
     m_TestVector->SetElement(2,c);
 
 
     output << "The vectorlength should be";
     output << sqrt(a*a+b*b+c*c);
     CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), sqrt(a*a+b*b+c*c), m_TestVector->GetNorm());
     output.flush();
 
     m_TestVector->Normalize();
 
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vectorlength should be 1.", 1.0, m_TestVector->GetNorm());
+    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vectorlength should be 1.", true, m_TestVector->GetNorm()-1 < DIF_VAL);
   }
 
   void TestRotatedVectorZeroDegrees()
   {
     int a = 1;
     int b = 2;
     int c = 3;
 
     double length;
 
     m_TestVector->SetElement(0,a);
     m_TestVector->SetElement(1,b);
     m_TestVector->SetElement(2,c);
 
     length = m_TestVector->GetNorm();
 
-    m_TestVector->Rotate(0,0,0);
+    m_TestVector->Rotate(0,0);
 
     CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", length, m_TestVector->GetNorm());
 
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be 1.0", 1.0, m_TestVector->GetElement(0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index1 should be 2.0", 2.0, m_TestVector->GetElement(1));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index2 should be 3.0", 3.0, m_TestVector->GetElement(2));
+    CPPUNIT_ASSERT_MESSAGE("The vector value at index0 should be 1.0", m_TestVector->GetElement(0) - 1 < DIF_VAL);
+    CPPUNIT_ASSERT_MESSAGE("The vector value at index1 should be 2.0", m_TestVector->GetElement(1) - 2 < DIF_VAL);
+    CPPUNIT_ASSERT_MESSAGE("The vector value at index2 should be 3.0", m_TestVector->GetElement(2) - 3 < DIF_VAL);
 
   }
 
   void TestRotatedVectorPositiveDegrees()
   {
-    DIF_VAL = 0.00001;
-    std::stringstream output;
-
-    int a = 1;
-    int b = 2;
-    int c = 3;
-
-    double length;
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    length = m_TestVector->GetNorm();
-
-    m_TestVector->Rotate(M_PI/4,0,0);
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << 1.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), 1.0, m_TestVector->GetElement(0));
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output << (-3.0/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(-3.0/sqrt(2.0)+sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << (3.0/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str() , true, abs(m_TestVector->GetElement(2)-(3.0/sqrt(2.0)+sqrt(2.0))) < DIF_VAL);
-    output.flush();
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    m_TestVector->Rotate(0,M_PI/4,0);
+    MITK_INFO << atan2(0, 0);
 
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << (2.0*sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(0)-(2.0*sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output <<2.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(2.0))          < DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << (2.0*sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(2)-(sqrt(2.0)))    < DIF_VAL);
-    output.flush();
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    m_TestVector->Rotate(0,0,M_PI/4);
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << (1/sqrt(2.0)-sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(0)-(1/sqrt(2.0)-sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output << (1/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(1/sqrt(2.0)+sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << 3.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(2)-(3.0))                  < DIF_VAL);
-    output.flush();
+    for(int r = 0; r<10; r++)
+    {
+      for(double phi = 0.1; phi < 3; phi+=0.1)
+      {
+        for(double theta = 0.1; theta < 3; theta+=0.1)
+        {
+
+          double rotateTheta = 0.1;
+          double rotatePhi = 0.1;
+
+          m_TestVector->SetElement(0, r * sin(theta) * cos(phi));
+          m_TestVector->SetElement(1, r * sin(theta) * sin(phi));
+          m_TestVector->SetElement(2, r * cos(theta));
+
+          m_TestVector->Rotate(rotateTheta, rotatePhi);
+
+          double newTheta = fmod(theta + rotateTheta, TWO_PI);
+          double newPhi = fmod(phi + rotatePhi, TWO_PI);
+
+          double expectedX = r * sin(newTheta) * cos(newPhi);
+          double expectedY = r * sin(newTheta) * sin(newPhi);
+          double expectedZ = r * cos(newTheta);
+
+          CPPUNIT_ASSERT_MESSAGE("The vector value at index0 should be " + std::to_string(expectedX) + " but was " + std::to_string(m_TestVector->GetElement(0))
+                                 + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta),
+                                 m_TestVector->GetElement(0) - expectedX < DIF_VAL);
+          CPPUNIT_ASSERT_MESSAGE("The vector value at index1 should be " + std::to_string(expectedY) + " but was " + std::to_string(m_TestVector->GetElement(0))
+                                 + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta),
+                                 m_TestVector->GetElement(1) - expectedY < DIF_VAL);
+          CPPUNIT_ASSERT_MESSAGE("The vector value at index2 should be " + std::to_string(expectedZ) + " but was " + std::to_string(m_TestVector->GetElement(0))
+                                 + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta),
+                                 m_TestVector->GetElement(2) - expectedZ < DIF_VAL);
+        }
+      }
+    }
 
   }
 
   void TestRotatedVectorNegativeDegrees()
   {
-    DIF_VAL = 0.00001;
-    std::stringstream output;
-
-    int a = 1;
-    int b = 2;
-    int c = 3;
 
-    double length;
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    length = m_TestVector->GetNorm();
-
-    m_TestVector->Rotate((-M_PI/4),0,0);
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << 1.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), 1.0,      m_TestVector->GetElement(0));
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output << (3.0/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(3.0/sqrt(2.0)+sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << (3.0/sqrt(2.0)-sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(2)-(3.0/sqrt(2.0)-sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    m_TestVector->Rotate(0,(-M_PI/4),0);
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << (-sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(0)-(-sqrt(2.0))) < DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output << 2.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(2.0)) < DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << (2.0 * sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(2)-(2.0 * sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    m_TestVector->Rotate(0,0,(-M_PI/4));
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
-
-    output << "The vector value at index0 should be";
-    output << (1/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(0)-(1/sqrt(2.0)+sqrt(2.0))) < DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index1 should be";
-    output << (-1/sqrt(2.0)+sqrt(2.0));
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(1)-(-1/sqrt(2.0)+sqrt(2.0)))< DIF_VAL);
-    output.flush();
-
-    output << "The vector value at index2 should be";
-    output << 3.0;
-    CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, abs(m_TestVector->GetElement(2)-(3.0)) < DIF_VAL);
-    output.flush();
+    //TODO rewrite this test
 
   }
 
   void TestRotatedVectorAroundAllAxis()
   {
-    DIF_VAL = 0.001;
-    int a = 1;
-    int b = 2;
-    int c = 3;
-
-    double length;
-
-    m_TestVector->SetElement(0,a);
-    m_TestVector->SetElement(1,b);
-    m_TestVector->SetElement(2,c);
-
-    length = m_TestVector->GetNorm();
-
-    m_TestVector->Rotate(M_PI/4, M_PI/4, M_PI/4);
-
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", true, abs(m_TestVector->GetNorm()-length) < DIF_VAL);
 
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be 2.76992", true, abs(m_TestVector->GetElement(0)-(2.76992))< DIF_VAL);
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index1 should be 1.77007", true, abs(m_TestVector->GetElement(1)-(1.77007))< DIF_VAL);
-    CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index2 should be 1.79605", true, abs(m_TestVector->GetElement(2)-(1.79605))< DIF_VAL);
+    //TODO rewrite this test
 
   }
 
   void TestScaledVector()
   {
     double a = 1.0;
     double b = 2.0;
     double c = 3.0;
 
     double length;
 
     for (double testFactor = -2.0; testFactor<=2.0; testFactor+=0.3)
     {
 
       double potElement0Fctr;
       double potElement1Fctr;
       double potElement2Fctr;
       double testLength;
 
       std::stringstream output;
 
       m_TestVector->SetElement(0,a);
       m_TestVector->SetElement(1,b);
       m_TestVector->SetElement(2,c);
 
       length = m_TestVector->GetNorm();
 
       potElement0Fctr = (m_TestVector->GetElement(0)*testFactor)*(m_TestVector->GetElement(0)*testFactor);
       potElement1Fctr = (m_TestVector->GetElement(1)*testFactor)*(m_TestVector->GetElement(1)*testFactor);
       potElement2Fctr = (m_TestVector->GetElement(2)*testFactor)*(m_TestVector->GetElement(2)*testFactor);
 
       testLength = sqrt(potElement0Fctr + potElement1Fctr + potElement2Fctr);
 
       m_TestVector->Scale(testFactor);
 
       CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should not be equal",
                                    sqrt(potElement0Fctr + potElement1Fctr + potElement2Fctr), m_TestVector->GetNorm());
 
       output << "The vector value at index0 should be";
       output << a*testFactor;
       CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), a*testFactor, m_TestVector->GetElement(0));
       output.flush();
 
       output << "The vector value at index1 should be";
       output << b*testFactor;
       CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), b*testFactor, m_TestVector->GetElement(1));
       output.flush();
 
       output << "The vector value at index2 should be";
       output << c*testFactor;
       CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), c*testFactor, m_TestVector->GetElement(2));
       output.flush();
     }
 
   }
 
   void TestClonedVector()
   {
     int a = 1;
     int b = 2;
     int c = 3;
 
     m_TestVector->SetElement(0,a);
     m_TestVector->SetElement(1,b);
     m_TestVector->SetElement(2,c);
 
     m_TestReturnVector = m_TestVector->Clone();
 
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector length should be equal",(m_TestVector->GetNorm()),m_TestReturnVector->GetNorm());
 
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be equal", m_TestVector->GetElement(0), m_TestReturnVector->GetElement(0));
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index1 should be equal", m_TestVector->GetElement(1), m_TestReturnVector->GetElement(1));
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index2 should be equal", m_TestVector->GetElement(2), m_TestReturnVector->GetElement(2));
 
-    m_TestReturnVector->Rotate(M_PI/4,M_PI/4,0);
+    m_TestReturnVector->Rotate(M_PI/4,M_PI/4);
 
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(0)!= m_TestReturnVector->GetElement(0));
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(1)!= m_TestReturnVector->GetElement(1));
     CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(2)!= m_TestReturnVector->GetElement(2));
 
     for(double testFactor = -2.0; testFactor<=2.0; testFactor+=0.3)
     {
 
       m_TestReturnVector->Scale(testFactor);
 
       CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(0)!= m_TestReturnVector->GetElement(0));
       CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(1)!= m_TestReturnVector->GetElement(1));
       CPPUNIT_ASSERT_EQUAL_MESSAGE( "The vector value at index0 should be not equal", true, m_TestVector->GetElement(2)!= m_TestReturnVector->GetElement(2));
     }
   }
 
   void tearDown()
   {
     m_TestVector = nullptr;
     m_TestReturnVector = nullptr;
   }
 
 };
 
 MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVector)
diff --git a/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVesselMeanderStrategyTest.cpp b/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVesselMeanderStrategyTest.cpp
index b427267f28..f6fb287d66 100644
--- a/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVesselMeanderStrategyTest.cpp
+++ b/Modules/PhotoacousticSimulation/Testing/mitkPhotoacousticVesselMeanderStrategyTest.cpp
@@ -1,127 +1,127 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkTestFixture.h>
 #include <mitkTestingMacros.h>
 
 #include <mitkImageReadAccessor.h>
 
 #include "mitkPhotoacousticVolume.h"
 #include "mitkPhotoacousticTissueGenerator.h"
 #include "mitkPhotoacousticVesselMeanderStrategy.h"
 
 class mitkPhotoacousticVesselMeanderStrategyTestSuite : public mitk::TestFixture
 {
   CPPUNIT_TEST_SUITE(mitkPhotoacousticVesselMeanderStrategyTestSuite);
   MITK_TEST(TestCalculateNewPositionInStraightLine);
   MITK_TEST(TestCalculateRandomlyDivergingPosition);
 
   CPPUNIT_TEST_SUITE_END();
 
 private:
   mitk::PhotoacousticVesselMeanderStrategy::Pointer m_TestVector;
   mitk::PhotoacousticSmartVector::Pointer m_TestPostion;
   mitk::PhotoacousticSmartVector::Pointer m_TestDirection;
 
 public:
 
   void setUp()
   {
 
     m_TestVector    = mitk::PhotoacousticVesselMeanderStrategy::New();
     m_TestPostion   = mitk::PhotoacousticSmartVector::New();
     m_TestDirection = mitk::PhotoacousticSmartVector::New();
 
   }
 
   void TestCalculateNewPositionInStraightLine()
   {
     std::stringstream output;
 
     int a = 0;
     int b = 1;
     int c = 2;
     int d = 3;
     int e = 4;
     int f = 5;
 
     for(int i = -2; i <= 2; i++)
     {
       if(i == 0)
       {
         i++;
       }
 
       for(int j = -2; j <= 2; j++ )
       {
         if(j == 0)
         {
           j++;
         }
 
         m_TestPostion -> SetElement(0,a*i);
         m_TestPostion -> SetElement(1,b*i);
         m_TestPostion -> SetElement(2,c*i);
         m_TestDirection -> SetElement(0,d*j);
         m_TestDirection -> SetElement(1,e*j);
         m_TestDirection -> SetElement(2,f*j);
 
-        m_TestVector->CalculateNewPositionInStraightLine(m_TestPostion,m_TestDirection,0, std::mt19937());
+        m_TestVector->CalculateNewPositionInStraightLine(m_TestPostion,m_TestDirection, 0, nullptr);
 
         MITK_INFO<<"m_TestPosition Element(0) ="<<m_TestPostion->GetElement(0);
         MITK_INFO<<"Data0 ="<<(a*i) + (d*j);
         MITK_INFO<<"m_TestPosition Element(1) ="<<m_TestPostion->GetElement(1);
         MITK_INFO<<"Data1 ="<<(b*i) + (e*j);
         MITK_INFO<<"m_TestPosition Element(2) ="<<m_TestPostion->GetElement(2);
         MITK_INFO<<"Data2 ="<<(c*i) + (f*j);
 
         output << "Element0 from m_TestPosition should be ";
         output << a*i + d*j;
         CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, a*i + d*j == m_TestPostion->GetElement(0));
         output.flush();
 
         output << "Element1 from m_TestPosition should be ";
         output << b*i + e*j;
         CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, b*i + e*j == m_TestPostion->GetElement(1));
         output.flush();
 
         output << "Element2 from m_TestPosition should be ";
         output << c*i + f*j;
         CPPUNIT_ASSERT_EQUAL_MESSAGE( output.str(), true, c*i + f*j == m_TestPostion->GetElement(2));
         output.flush();
 
       }
     }
 
   }
 
   void TestCalculateRandomlyDivergingPosition()
   {
 
   }
 
   void tearDown()
   {
 
     m_TestVector    = nullptr;
     m_TestPostion   = nullptr;
     m_TestDirection = nullptr;
 
   }
 
 };
 
 MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVesselMeanderStrategy)
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp
index 9ada0d316b..6a1df077c2 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp
@@ -1,63 +1,45 @@
 #include "mitkPhotoacousticFatTissueGenerator.h"
 #include <random>
 #include <chrono>
 #include <mitkCommon.h>
 
 void mitk::PhotoacousticFatTissueGenerator::CreateFatGlobules(mitk::PhotoacousticVolume::Pointer image, double bodyFatPercentage,
                                                               double fatAbsorption, double fatScattering, double fatAnisotropy,
-                                                              bool useRngSeed, long rngSeed)
+                                                              std::mt19937* rng)
 {
-  std::mt19937 rng;
-  std::random_device randomDevice;
-  if(useRngSeed)
-  {
-    rng.seed(rngSeed);
-  }
-  else
-  {
-    if(randomDevice.entropy()>0.1)
-    {
-      rng.seed(randomDevice());
-    }
-    else
-    {
-      rng.seed(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
-    }
-  }
-
   int xDim = image->getXDim();
   int yDim = image->getYDim();
   int zDim = image->getZDim();
 
   long length = xDim*yDim*zDim;
 
-  double sqrtLength = sqrt(sqrt(length));
+  double sqrtLength = 2 * sqrt(sqrt(length));
 
   std::uniform_int_distribution<> numberOfFatGlobulesDist(0.83*sqrtLength, 1.2*sqrtLength);
-  std::normal_distribution<> normal(1, 0.05);
+  std::normal_distribution<> normal(1, 0.1);
 
-  int numberOfFatGlobules = numberOfFatGlobulesDist(rng);
+  int numberOfFatGlobules = numberOfFatGlobulesDist(*rng);
 
   double radius = pow(length*bodyFatPercentage/100.0/(numberOfFatGlobules * 4.0 / 3.0 * 3.141), (1.0/3.0));
 
   std::uniform_int_distribution<> xDist(-radius, xDim+radius);
   std::uniform_int_distribution<> yDist(-radius, yDim+radius);
   std::uniform_int_distribution<> zDist(1.5*radius, zDim+radius);
 
   for(int fatGlobuleIdx = 0; fatGlobuleIdx<numberOfFatGlobules; ++fatGlobuleIdx)
   {
-    double usedradius = radius * normal(rng);
-    int globX = xDist(rng);
-    int globY = yDist(rng);
-    int globZ = zDist(rng);
+    double usedradius = radius * normal(*rng);
+    int globX = xDist(*rng);
+    int globY = yDist(*rng);
+    int globZ = zDist(*rng);
     for(int x = 0; x < xDim; x++)
       for(int y = 0; y < yDim; y++)
         for(int z = 0; z < zDim; z++)
         {
           if((globX - x) * (globX - x) + (globY - y) * (globY - y) + (globZ - z) * (globZ - z) <= usedradius * usedradius)
           {
             image->SetVolumeValues(x, y, z, fatAbsorption, fatScattering, fatAnisotropy, mitk::PhotoacousticVolume::SegmentationType::FAT);
           }
         }
   }
 }
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h
index 9dde4dabc0..6509b43454 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h
@@ -1,41 +1,41 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef MITKPHOTOACOUSTICFATTISSUEGENERATOR_H
 #define MITKPHOTOACOUSTICFATTISSUEGENERATOR_H
 
 #include "MitkPhotoacousticSimulationExports.h"
 
 #include "mitkPhotoacousticVolume.h"
 
 namespace mitk
 {
 class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticFatTissueGenerator
 {
 public:
 
   static void CreateFatGlobules(mitk::PhotoacousticVolume::Pointer image, double bodyFatPercentage,
                                double fatAbsorption, double fatScattering, double fatAnisotropy,
-                                bool useRngSeed, long rngSeed);
+                                std::mt19937* rng);
 
 private:
   PhotoacousticFatTissueGenerator();
   ~PhotoacousticFatTissueGenerator();
 
 };
 }
 
 #endif // MITKPHOTOACOUSTICFATTISSUEGENERATOR_H
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp
index dfa6d2b2ff..5e310f088f 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp
@@ -1,163 +1,162 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticTissueGenerator.h"
 #include "mitkPhotoacousticSmartVector.h"
 #include "mitkPhotoacousticFatTissueGenerator.h"
 
 mitk::PhotoacousticVolume::Pointer mitk::PhotoacousticTissueGenerator::GenerateInSilicoData(int xDim, int yDim, int zDim, double spacing,
                                                                                             double basicAbsorption, double basicScattering, double basicAnisotropy, double airThickness,
                                                                                             double skinAbsorption, double skinScattering, double skinAnisotropy, double skinThickness,
                                                                                             double fatAbsorption, double fatScattering, double fatAnisotropy, double fatPercentage,
                                                                                             bool randomizePhysicalProperties, double randomPercentage,
                                                                                             int minNumberOfVessels, int maxNumberOfVessels,
                                                                                             double minBending, double maxBending,
                                                                                             double minAbs, double maxAbs,
                                                                                             double minRadius, double maxRadius,
                                                                                             double mcflag, double launchflag, double boundaryflag,
                                                                                             double launchPointX, double launchPointY, double launchPointZ,
                                                                                             double focusPointX, double focusPointY, double focusPointZ,
                                                                                             double trajectoryVectorX, double trajectoryVectorY, double trajectoryVectorZ,
                                                                                             double radius, double waist,
                                                                                             double sigma, bool doGauss, int bifurcationFrequency,
                                                                                             double vesselScattering, double vesselAnisotropy,
                                                                                             bool useRngSeed, long rngSeed,
                                                                                             mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewVesselPositionCallback)
 {
   MITK_DEBUG << "Initializing Empty Volume";
   mitk::PhotoacousticVolume::Pointer generatedVolume = mitk::PhotoacousticVolume::New();
   generatedVolume->Initialize(xDim, yDim, zDim, spacing,
                               basicAbsorption, basicScattering, basicAnisotropy,
                               mcflag, launchflag, boundaryflag,
                               launchPointX, launchPointY, launchPointZ,
                               focusPointX, focusPointY, focusPointZ,
                               trajectoryVectorX, trajectoryVectorY, trajectoryVectorZ,
                               radius, waist,
                               sigma, doGauss,
                               airThickness, skinThickness,
                               skinAbsorption, skinScattering, skinAnisotropy);
 
-  mitk::PhotoacousticFatTissueGenerator::CreateFatGlobules(generatedVolume, fatPercentage, fatAbsorption, fatScattering, fatAnisotropy, useRngSeed, rngSeed);
-
   std::mt19937 randomNumberGenerator;
   std::random_device randomDevice;
   if(useRngSeed)
   {
     randomNumberGenerator.seed(rngSeed);
   }
   else
   {
     if(randomDevice.entropy()>0.1)
     {
       randomNumberGenerator.seed(randomDevice());
     }
     else
     {
       randomNumberGenerator.seed(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
     }
   }
 
+  mitk::PhotoacousticFatTissueGenerator::CreateFatGlobules(generatedVolume, fatPercentage, fatAbsorption, fatScattering, fatAnisotropy, &randomNumberGenerator);
+
   std::uniform_int_distribution<int> randomNumVesselDistribution(minNumberOfVessels, maxNumberOfVessels);
   std::uniform_real_distribution<double> randomBendingDistribution(minBending, maxBending);
   std::uniform_real_distribution<double> randomAbsorptionDistribution(minAbs, maxAbs);
   std::uniform_real_distribution<double> randomRadiusDistribution(minRadius, maxRadius);
   std::uniform_real_distribution<double> randomSign(-1, 1);
 
   int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator);
 
   generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels);
   generatedVolume->AddIntProperty("bifurcationFrequency", bifurcationFrequency);
 
   MITK_DEBUG << "Simulating vessel structures";
   for(int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++)
   {
 
     mitk::PhotoacousticSmartVector::Pointer initialPosition = mitk::PhotoacousticSmartVector::New();
     mitk::PhotoacousticSmartVector::Pointer initialDirection = mitk::PhotoacousticSmartVector::New();
 
     double initialRadius = randomRadiusDistribution(randomNumberGenerator) / spacing / 10;
     std::stringstream radiusString;
     radiusString <<"vessel_"<<vesselNumber+1<<"_radius";
     generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius);
 
     double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator);
     std::stringstream absorptionString;
     absorptionString <<"vessel_"<<vesselNumber+1<<"_absorption";
     generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient);
 
     double bendingFactor = randomBendingDistribution(randomNumberGenerator);
     std::stringstream bendingString;
     bendingString <<"vessel_"<<vesselNumber+1<<"_bendingFactor";
     generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor);
 
     bool stickXBorder = randomSign(randomNumberGenerator) < 0;
     bool stickLeft = randomSign(randomNumberGenerator) < 0;
 
     if(stickXBorder)
     {
       if(stickLeft)
       {
-        initialPosition->Randomize(0, 0, 0, yDim, 3+initialRadius, zDim, randomNumberGenerator);
-        initialDirection->Randomize(1, 2, -0.2, 0.2, -0.2, 0.2,randomNumberGenerator);
+        initialPosition->Randomize(0, 0, 0, yDim, 3+initialRadius, zDim, &randomNumberGenerator);
+        initialDirection->Randomize(1, 2, -0.2, 0.2, -0.2, 0.2,&randomNumberGenerator);
       }
       else
       {
-        initialPosition->Randomize(xDim, xDim, 0, yDim, 3+initialRadius, zDim, randomNumberGenerator);
-        initialDirection->Randomize(-2, -1, -0.2, 0.2, -0.2, 0.2, randomNumberGenerator);
+        initialPosition->Randomize(xDim, xDim, 0, yDim, 3+initialRadius, zDim, &randomNumberGenerator);
+        initialDirection->Randomize(-2, -1, -0.2, 0.2, -0.2, 0.2, &randomNumberGenerator);
       }
     }
     else
     {
       if(stickLeft)
       {
-        initialPosition->Randomize(0, xDim, 0, 0, 3+initialRadius, zDim, randomNumberGenerator);
-        initialDirection->Randomize(-0.2, 0.2, 1, 2, -0.2, 0.2, randomNumberGenerator);
+        initialPosition->Randomize(0, xDim, 0, 0, 3+initialRadius, zDim, &randomNumberGenerator);
+        initialDirection->Randomize(-0.2, 0.2, 1, 2, -0.2, 0.2, &randomNumberGenerator);
       }
       else
       {
-        initialPosition->Randomize(0, xDim, yDim, yDim, 3+initialRadius, zDim, randomNumberGenerator);
-        initialDirection->Randomize(-0.2, 0.2, -2, -1, -0.2, 0.2, randomNumberGenerator);
+        initialPosition->Randomize(0, xDim, yDim, yDim, 3+initialRadius, zDim, &randomNumberGenerator);
+        initialDirection->Randomize(-0.2, 0.2, -2, -1, -0.2, 0.2, &randomNumberGenerator);
       }
     }
 
     mitk::PhotoacousticVesselTree::Pointer vesselStructure = mitk::PhotoacousticVesselTree::New();
-    vesselStructure->SetBifurcationFrequency(bifurcationFrequency);
 
     std::stringstream scatteringString;
     scatteringString <<"vessel_"<<vesselNumber+1<<"_scattering";
     generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering);
 
     std::stringstream anisotropyString;
     anisotropyString <<"vessel_"<<vesselNumber+1<<"_anisotropy";
     generatedVolume->AddDoubleProperty(anisotropyString.str(), vesselAnisotropy);
 
     MITK_DEBUG << "Initializing vessel structure";
-    vesselStructure->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, vesselScattering, vesselAnisotropy);
+    vesselStructure->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, vesselScattering, vesselAnisotropy, bifurcationFrequency, &randomNumberGenerator);
 
     MITK_DEBUG << "Simulating vessel structure";
     while(!vesselStructure->IsFinished())
     {
-      vesselStructure->Step(generatedVolume, calculateNewVesselPositionCallback, bendingFactor, randomNumberGenerator);
+      vesselStructure->Step(generatedVolume, calculateNewVesselPositionCallback, bendingFactor, &randomNumberGenerator);
     }
 
   }
 
   generatedVolume->GaussianBlur3D();
 
   generatedVolume->FinalizeVolume(randomizePhysicalProperties, randomPercentage, useRngSeed, rngSeed);
 
   return generatedVolume;
 }
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp
index 0db95091ff..ca4c4efb5c 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp
@@ -1,156 +1,162 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticVessel.h"
 #include <chrono>
 #include <math.h>
 #include <thread>
 
-mitk::PhotoacousticVessel::PhotoacousticVessel()
+mitk::PhotoacousticVessel::PhotoacousticVessel():
+  m_Range(M_PI/16, M_PI/8),
+  m_Sign(-1, 1)
 {
   srand(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
   m_Finished = false;
   m_Position = nullptr;
   m_Direction = nullptr;
   m_Radius = -1;
   m_AbsorptionCoefficient = -1;
   m_ScatteringCoefficient = -1;
   m_AnisotropyCoefficient = -1;
   m_VesselMeanderStrategy = nullptr;
   m_WalkedDistance = 0;
   m_State = UNINITIALIZED;
 }
 
 mitk::PhotoacousticVessel::~PhotoacousticVessel()
 {
   m_Position = nullptr;
   m_Direction = nullptr;
   m_VesselMeanderStrategy = nullptr;
   m_State = UNINITIALIZED;
 }
 
 void mitk::PhotoacousticVessel::Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition,
                                            mitk::PhotoacousticSmartVector::Pointer initialDirection,
                                            double initialRadius, double absorptionCoefficient,
-                                           double scatteringCoefficient, double anisotophyCoefficient)
+                                           double scatteringCoefficient, double anisotophyCoefficient,
+                                           double bifurcationFrequency, std::mt19937* rng)
 {
-  m_Position = initialPosition;
-  m_Direction = initialDirection;
   m_Radius = initialRadius;
   m_AbsorptionCoefficient = absorptionCoefficient;
   m_ScatteringCoefficient = scatteringCoefficient;
   m_AnisotropyCoefficient = anisotophyCoefficient;
-  m_VesselMeanderStrategy = mitk::PhotoacousticVesselMeanderStrategy::New();
   m_Finished = false;
+  m_Position = initialPosition;
+  m_Direction = initialDirection;
+  m_VesselMeanderStrategy = mitk::PhotoacousticVesselMeanderStrategy::New();
   m_State = INITIALIZED;
+  std::normal_distribution<> factor(1, 0.2);
+  m_BifurcationFrequency = bifurcationFrequency * factor(*rng);
 }
 
-void mitk::PhotoacousticVessel::ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937 rng)
+void mitk::PhotoacousticVessel::ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng)
 {
   AssertState(INITIALIZED);
   mitk::PhotoacousticSmartVector::Pointer oldPosition = mitk::PhotoacousticSmartVector::New();
   oldPosition->SetValue(m_Position);
 
   (m_VesselMeanderStrategy->*calculateNewPosition)(m_Position, m_Direction, bendingFactor, rng);
 
   DrawVesselInVolume(oldPosition, volume);
 }
 
-bool mitk::PhotoacousticVessel::CanBifurcate(double bifurcationFrequency)
+bool mitk::PhotoacousticVessel::CanBifurcate()
 {
-  return bifurcationFrequency < m_WalkedDistance;
+  return m_BifurcationFrequency < m_WalkedDistance;
 }
 
-mitk::PhotoacousticVessel::Pointer mitk::PhotoacousticVessel::Bifurcate(std::mt19937 rng)
+mitk::PhotoacousticVessel::Pointer mitk::PhotoacousticVessel::Bifurcate(std::mt19937* rng)
 {
   AssertState(INITIALIZED);
   mitk::PhotoacousticVessel::Pointer newVessel = mitk::PhotoacousticVessel::New();
   mitk::PhotoacousticSmartVector::Pointer newDirection = m_Direction->Clone();
 
-  std::uniform_real_distribution<> range(M_PI/16, M_PI/8);
+  double thetaChange = m_Range(*rng) * (m_Sign(*rng) < 1 ? -1 : 1);
+  double phiChange = m_Range(*rng) * (m_Sign(*rng) < 1 ? -1 : 1);
 
-  double xAngle = range(rng);
-  double yAngle = range(rng);
-  double zAngle = range(rng);
+  MITK_INFO << "Theta change: " << thetaChange;
+  MITK_INFO << "Phi change: " << phiChange;
+  MITK_INFO << "(m_Sign(*rng) = " << (m_Sign(*rng);
 
-  newDirection->Rotate(xAngle, yAngle, zAngle);
+  newDirection->Rotate(thetaChange, phiChange);
 
-  m_Direction->Rotate(-xAngle, -yAngle, -zAngle);
+  m_Direction->Rotate(-thetaChange, -phiChange);
 
   std::uniform_real_distribution<> rangeRadius(m_Radius*0.6, m_Radius*0.8);
 
-  double newRadius = rangeRadius(rng);
+  double newRadius = rangeRadius(*rng);
 
   m_Radius = sqrt(m_Radius*m_Radius - newRadius*newRadius);
 
-  newVessel->Initialize(m_Position->Clone(), newDirection, newRadius, m_AbsorptionCoefficient, m_ScatteringCoefficient, m_AnisotropyCoefficient);
+  newVessel->Initialize(m_Position->Clone(), newDirection, newRadius, m_AbsorptionCoefficient, m_ScatteringCoefficient, m_AnisotropyCoefficient, m_BifurcationFrequency, rng);
 
   m_WalkedDistance = 0;
 
   return newVessel;
 }
 
 void mitk::PhotoacousticVessel::DrawVesselInVolume(mitk::PhotoacousticSmartVector::Pointer fromPosition, mitk::PhotoacousticVolume::Pointer volume)
 {
   AssertState(INITIALIZED);
   mitk::PhotoacousticSmartVector::Pointer diffVector = mitk::PhotoacousticSmartVector::New();
   diffVector->SetElement(0, fromPosition->GetElement(0)-m_Position->GetElement(0));
   diffVector->SetElement(1, fromPosition->GetElement(1)-m_Position->GetElement(1));
   diffVector->SetElement(2, fromPosition->GetElement(2)-m_Position->GetElement(2));
 
   mitk::PhotoacousticSmartVector::Pointer stepSize = mitk::PhotoacousticSmartVector::New();
   stepSize->SetValue(m_Direction);
   stepSize->Scale(SCALING_FACTOR);
 
   while (diffVector->GetNorm() >= SCALING_FACTOR)
   {
     m_WalkedDistance += stepSize->GetNorm();
 
     fromPosition->SetElement(0, fromPosition->GetElement(0)+stepSize->GetElement(0));
     fromPosition->SetElement(1, fromPosition->GetElement(1)+stepSize->GetElement(1));
     fromPosition->SetElement(2, fromPosition->GetElement(2)+stepSize->GetElement(2));
 
     int xPos = fromPosition->GetElement(0);
     int yPos = fromPosition->GetElement(1);
     int zPos = fromPosition->GetElement(2);
 
     for(int x = xPos - m_Radius; x <= xPos + m_Radius; x++)
       for(int y = yPos - m_Radius; y <= yPos + m_Radius; y++)
         for(int z = zPos - m_Radius; z <= zPos + m_Radius; z++)
         {
           if( m_Radius*m_Radius >= (x-xPos)*(x-xPos)+(y-yPos)*(y-yPos)+(z-zPos)*(z-zPos))
           {
             volume->SetVolumeValues(x, y, z, m_AbsorptionCoefficient, m_ScatteringCoefficient, m_AnisotropyCoefficient, mitk::PhotoacousticVolume::SegmentationType::VESSEL);
           }
         }
 
     diffVector->SetElement(0, fromPosition->GetElement(0)-m_Position->GetElement(0));
     diffVector->SetElement(1, fromPosition->GetElement(1)-m_Position->GetElement(1));
     diffVector->SetElement(2, fromPosition->GetElement(2)-m_Position->GetElement(2));
 
     if(!volume->IsInsideVolume(xPos, yPos, zPos))
     {
       m_Radius = 0;
       break;
     }
   }
 
 }
 
 bool mitk::PhotoacousticVessel::IsFinished()
 {
   return m_Radius < MIN_RADIUS;
 }
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h
index 5b5c9ac7d4..0e6364acfd 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h
@@ -1,113 +1,119 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef MITKVESSEL_H
 #define MITKVESSEL_H
 
 #include "mitkVector.h"
 #include "mitkPhotoacousticVesselMeanderStrategy.h"
 #include "mitkPhotoacousticVolume.h"
 #include "mitkPhotoacousticSmartVector.h"
 #include "mitkPhotoacousticStatefulObject.h"
 
 #include <MitkPhotoacousticSimulationExports.h>
 
 //Includes for smart pointer usage
 #include "mitkCommon.h"
 #include "itkLightObject.h"
 
 namespace mitk {
 
 class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticVessel : public itk::LightObject, public mitk::PhotoacousticStatefulObject
 {
 public:
   mitkClassMacroItkParent(mitk::PhotoacousticVessel, itk::LightObject)
   itkFactorylessNewMacro(Self)
 
   /**
    * Callback function definition of a PhotoacousticVesselMeanderStrategy
    */
   typedef void (mitk::PhotoacousticVesselMeanderStrategy::*CalculateNewVesselPositionCallback)
-  (mitk::PhotoacousticSmartVector::Pointer, mitk::PhotoacousticSmartVector::Pointer, double, std::mt19937);
+  (mitk::PhotoacousticSmartVector::Pointer, mitk::PhotoacousticSmartVector::Pointer, double, std::mt19937*);
 
   /**
    * @brief Initialize initializes this PhotoacousticVessel with the given parameters.
    *
    * @param initialPosition
    * @param initialDirection
    * @param initialRadius
    * @param absorptionCoefficient
    * @param scatteringCoefficient
    * @param anisotropyCoefficient
+   * @param bifurcationFrequency
+   * @param rng
    */
   void Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition, mitk::PhotoacousticSmartVector::Pointer initialDirection,
-                  double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotropyCoefficient);
+                  double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotropyCoefficient,
+                  double bifurcationFrequency, std::mt19937* rng);
 
   /**
    * @brief ExpandVessel makes this PhotoacousticVessel expand one step in its current direction.
    * After expanding, the vessel will draw itself into the given PhotoacousticVolume.
    *
    * @param volume volume for the vessel to draw itself in
    * @param calculateNewPosition a callback function of the PhotoacousticVesselMeanderStrategy class.
    * It is used to  calculate the final position after taking the step.
    * @param bendingFactor a metric of how much the PhotoacousticVessel should bend. If set to 0 the vessel will go in a straight line.
    */
-  void ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937 rng);
+  void ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng);
 
   /**
    * @brief CanBifurcate
-   * @param bifurcationProbability
    * @return true if the PhotoacousticVessel is ready to Bifurcate()
    */
-  bool CanBifurcate(double bifurcationProbability);
+  bool CanBifurcate();
 
   /**
    * @brief Bifurcate bifurcates this vessel into two new ones. Makes sure that the volume of the vessels stays the same.
    *
    * @return a new vessel split up from the current one.
    */
-  mitk::PhotoacousticVessel::Pointer Bifurcate(std::mt19937 rng);
+  mitk::PhotoacousticVessel::Pointer Bifurcate(std::mt19937* rng);
 
   /**
    * @brief IsFinished
    * @return true if the vessel cannot expand any further
    */
   bool IsFinished();
 
 protected:
   PhotoacousticVessel();
   virtual ~PhotoacousticVessel();
 
 private:
 
   const double MIN_RADIUS = 1;
   const double SCALING_FACTOR = 0.33;
 
   void DrawVesselInVolume(mitk::PhotoacousticSmartVector::Pointer toPosition, mitk::PhotoacousticVolume::Pointer volume);
   mitk::PhotoacousticSmartVector::Pointer m_Position;
   mitk::PhotoacousticSmartVector::Pointer m_Direction;
   double m_Radius;
   double m_AbsorptionCoefficient;
   double m_ScatteringCoefficient;
   double m_AnisotropyCoefficient;
   double m_WalkedDistance;
+  double m_BifurcationFrequency;
   mitk::PhotoacousticVesselMeanderStrategy::Pointer m_VesselMeanderStrategy;
   bool m_Finished;
+
+  std::uniform_real_distribution<> m_Range;
+  std::uniform_real_distribution<> m_Sign;
 };
 
 }
 
 #endif // MITKVESSEL_H
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp
index 75886ea4a4..63a4e803a5 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp
@@ -1,53 +1,53 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticVesselMeanderStrategy.h"
 
 mitk::PhotoacousticVesselMeanderStrategy::PhotoacousticVesselMeanderStrategy()
 {
 }
 
 mitk::PhotoacousticVesselMeanderStrategy::~PhotoacousticVesselMeanderStrategy()
 {
 }
 
-void mitk::PhotoacousticVesselMeanderStrategy::CalculateNewPositionInStraightLine(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double /*bendingFactor*/, std::mt19937 rng)
+void mitk::PhotoacousticVesselMeanderStrategy::CalculateNewPositionInStraightLine(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double /*bendingFactor*/, std::mt19937* rng)
 {
   if(direction->GetNorm() <= mitk::eps)
   {
     direction->Randomize(rng);
   }
 
   position->SetElement(0,position->GetElement(0) + direction->GetElement(0));
   position->SetElement(1,position->GetElement(1) + direction->GetElement(1));
   position->SetElement(2,position->GetElement(2) + direction->GetElement(2));
 }
 
-void mitk::PhotoacousticVesselMeanderStrategy::CalculateRandomlyDivergingPosition(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937 rng)
+void mitk::PhotoacousticVesselMeanderStrategy::CalculateRandomlyDivergingPosition(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937* rng)
 {
   if(direction->GetNorm() <= mitk::eps)
   {
     direction->Randomize(rng);
   }
 
   direction->RandomizeByPercentage(0.1, bendingFactor, rng);
 
   position->SetElement(0,position->GetElement(0) + direction->GetElement(0));
   position->SetElement(1,position->GetElement(1) + direction->GetElement(1));
   position->SetElement(2,position->GetElement(2) + direction->GetElement(2));
 }
 
 
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h
index 3ff49e03aa..e78ef6f81c 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h
@@ -1,62 +1,62 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef MITKVESSELMEANDERSTRATEGY_H
 #define MITKVESSELMEANDERSTRATEGY_H
 
 #include "mitkVector.h"
 #include "mitkPhotoacousticSmartVector.h"
 
 #include <MitkPhotoacousticSimulationExports.h>
 
 //Includes for smart pointer usage
 #include "mitkCommon.h"
 #include "itkLightObject.h"
 
 namespace mitk {
 
 class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticVesselMeanderStrategy : public itk::LightObject
 {
 public:
   mitkClassMacroItkParent(mitk::PhotoacousticVesselMeanderStrategy, itk::LightObject)
   itkFactorylessNewMacro(Self)
 
   /**
    * @brief CalculateNewPositionInStraightLine calculates the new position by just following the direction vector.
    * @param position
    * @param direction
    * @param bendingFactor
    */
-  void CalculateNewPositionInStraightLine(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937 rng);
+  void CalculateNewPositionInStraightLine(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937* rng);
 
   /**
    * @brief CalculateRandomlyDivergingPosition calculates the new position by modifying the direction vector randomly,
    * proportional to the selected bendingFactor. This means, that the vessels will bend in each expansion step,
    * if bendingFactor > 0.
    *
    * @param position
    * @param direction
    * @param bendingFactor
    */
-  void CalculateRandomlyDivergingPosition(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937 rng);
+  void CalculateRandomlyDivergingPosition(mitk::PhotoacousticSmartVector::Pointer position, mitk::PhotoacousticSmartVector::Pointer direction, double bendingFactor, std::mt19937* rng);
 
 protected:
   PhotoacousticVesselMeanderStrategy();
   virtual ~PhotoacousticVesselMeanderStrategy();
 };
 }
 
 #endif // MITKVESSELMEANDERSTRATEGY_H
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp
index 040847a134..acb94954da 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp
@@ -1,84 +1,76 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticVesselTree.h"
 #include "thread"
 void mitk::PhotoacousticVesselTree::Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition, mitk::PhotoacousticSmartVector::Pointer initialDirection,
-                                               double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotopyCoefficient)
+                                               double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotopyCoefficient,
+                                               double bifurcationFrequency, std::mt19937* rng)
 {
   m_CurrentSubvessels = new std::vector<mitk::PhotoacousticVessel::Pointer>();
   mitk::PhotoacousticVessel::Pointer tmpVessel = mitk::PhotoacousticVessel::New();
-  tmpVessel->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, scatteringCoefficient, anisotopyCoefficient);
+  tmpVessel->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, scatteringCoefficient, anisotopyCoefficient,
+                        bifurcationFrequency, rng);
   m_CurrentSubvessels->push_back(tmpVessel);
   m_State = INITIALIZED;
 }
 
 mitk::PhotoacousticVesselTree::PhotoacousticVesselTree()
 {
   m_CurrentSubvessels = nullptr;
-  m_TmpCounter = 0;
-  m_BifurcationFrequency = 0.1;
   m_State = UNINITIALIZED;
 }
 
 mitk::PhotoacousticVesselTree::~PhotoacousticVesselTree()
 {
   m_CurrentSubvessels = nullptr;
-  m_TmpCounter = 0;
-  m_BifurcationFrequency = 0;
-  m_State = UNINITIALIZED;
 }
 
-void mitk::PhotoacousticVesselTree::Step(mitk::PhotoacousticVolume::Pointer volume, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937 rng)
+void mitk::PhotoacousticVesselTree::Step(mitk::PhotoacousticVolume::Pointer volume, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng)
 {
   AssertState(INITIALIZED);
   std::vector<mitk::PhotoacousticVessel::Pointer> newVessels;
 
   for(unsigned int vesselTreeIndex = 0; vesselTreeIndex < m_CurrentSubvessels->size(); vesselTreeIndex++)
   {
     mitk::PhotoacousticVessel::Pointer currentVessel = m_CurrentSubvessels->at(vesselTreeIndex);
     if(!currentVessel->IsFinished())
     {
       currentVessel->ExpandVessel(volume, calculateNewPosition, bendingFactor, rng);
-      if(currentVessel->CanBifurcate(m_BifurcationFrequency))
+      if(currentVessel->CanBifurcate())
       {
         newVessels.push_back(currentVessel->Bifurcate(rng));
       }
     }
   }
 
   for(unsigned int newVesselsIndex = 0; newVesselsIndex < newVessels.size(); newVesselsIndex++)
   {
     m_CurrentSubvessels->push_back(newVessels.at(newVesselsIndex));
   }
 
 }
 
 bool mitk::PhotoacousticVesselTree::IsFinished()
 {
   AssertState(INITIALIZED);
   for(unsigned int vesselTreeIndex = 0; vesselTreeIndex < m_CurrentSubvessels->size(); vesselTreeIndex++)
   {
     if(!m_CurrentSubvessels->at(vesselTreeIndex)->IsFinished())
       return false;
   }
   return true;
 }
-
-void mitk::PhotoacousticVesselTree::SetBifurcationFrequency(double frequency)
-{
-  m_BifurcationFrequency = frequency;
-}
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h
index 48350438f2..88a591fbbc 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h
@@ -1,89 +1,83 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef MITKVESSELSTRUCTURE_H
 #define MITKVESSELSTRUCTURE_H
 
 //std includes
 #include <vector>
 
 //mitk includes
 #include "mitkVector.h"
 #include "mitkPhotoacousticVessel.h"
 #include "mitkPhotoacousticVolume.h"
 #include "mitkPhotoacousticStatefulObject.h"
 #include <MitkPhotoacousticSimulationExports.h>
 
 //Includes for smart pointer usage
 #include "mitkCommon.h"
 #include "itkLightObject.h"
 
 namespace mitk {
 
 class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticVesselTree : public itk::LightObject, public mitk::PhotoacousticStatefulObject
 {
 public:
   mitkClassMacroItkParent(mitk::PhotoacousticVesselTree, itk::LightObject)
   itkFactorylessNewMacro(Self)
 
   /**
    * @brief Initializes this PhotoacousticVesselTree with the given parameters.
    *
    * @param initialPosition
    * @param initialDirection
    * @param initialRadius
    * @param absorbtionCoefficient
    * @param scatteringCoefficient
    * @param anisotropyCoefficient
+   * @param bifurcationFrequency
+   * @param rng
    */
   void Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition, mitk::PhotoacousticSmartVector::Pointer initialDirection,
-                  double initialRadius, double absorbtionCoefficient, double scatteringCoefficient, double anisotropyCoefficient);
+                  double initialRadius, double absorbtionCoefficient, double scatteringCoefficient, double anisotropyCoefficient,
+                  double bifurcationFrequency, std::mt19937* rng);
 
   /**
    * @brief Step Performs a simulation step, in which all subvessels of this PhotoacousticVesselTree are expanded.
    *
    * @param volume
    * @param calculateNewPosition
    * @param bendingFactor
    */
-  void Step(mitk::PhotoacousticVolume::Pointer volume, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937 rng);
+  void Step(mitk::PhotoacousticVolume::Pointer volume, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng);
 
   /**
    * @brief IsFinished
    * @return true if no subvessel can be expanded.
    */
   bool IsFinished();
 
-  /**
-   * @brief SetBifurcationFrequency
-   * @param frequency
-   */
-  void SetBifurcationFrequency(double frequency);
-
 protected:
   PhotoacousticVesselTree();
   virtual ~PhotoacousticVesselTree();
 
 private:
   std::vector<mitk::PhotoacousticVessel::Pointer>* m_CurrentSubvessels;
-  int m_TmpCounter;
-  double m_BifurcationFrequency;
-
 };
 
 }
 
 #endif // MITKVESSELSTRUCTURE_H
diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp
index c42a3e32ab..8880bef751 100644
--- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp
+++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp
@@ -1,451 +1,451 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include <mitkPhotoacousticVolume.h>
 #include <mitkProperties.h>
 #include <mitkCoreServices.h>
 #include <mitkIPropertyPersistence.h>
 #include <itkImage.h>
 #include <itkDiscreteGaussianImageFilter.h>
 #include <mitkImageCast.h>
 #include <mitkImageToItk.h>
 #include <chrono>
 
 void mitk::PhotoacousticVolume::Initialize(int xDim, int yDim, int zDim, double spacing, double initAbsorption, double initScattering, double initAnisotropy,
                                            int mcflag, int launchflag, int boundaryflag,
                                            double launchPointX, double launchPointY, double launchPointZ,
                                            double focusPointX, double focusPointY, double focusPointZ,
                                            double trajectoryVectorX, double trajectoryVectorY, double trajectoryVectorZ,
                                            double radius, double waist,
                                            double sigma, bool doGauss,
                                            double airThickness, double epidermisThickness,
                                            double skinAbsorption, double skinScattering, double skinAnisotropy)
 {
   AssertState(UNINITIALIZED);
 
   int size = xDim*yDim*zDim;
   m_AbsorptionArray = new double[size];
   m_ScatteringArray = new double[size];
   m_AnisotropyArray = new double[size];
   m_SegmentationArray = new double[size];
 
   m_ThicknessAir = airThickness;
   m_ThicknessSkin = epidermisThickness;
 
   m_AbsorptionSkin = skinAbsorption;
   m_ScatteringSkin = skinScattering;
   m_AnisotropySkin = skinAnisotropy;
 
 
   m_XDim = xDim;
   m_YDim = yDim;
   m_ZDim = zDim;
   m_TDim = 4;
 
   m_Spacing = spacing;
 
   for (int index = 0; index < size; index++)
   {
     m_AbsorptionArray[index]=initAbsorption;
     m_ScatteringArray[index]=initScattering;
     m_AnisotropyArray[index]=initAnisotropy;
     m_SegmentationArray[index] = SegmentationType::BACKGROUND;
   }
 
   if(doGauss)
   {
     m_Sigma = sigma;
   }
   else
   {
     m_Sigma = 0;
   }
 
   //Set properties
   AddIntProperty("mcflag", mcflag);
   AddIntProperty("launchflag", launchflag);
   AddIntProperty("boundaryflag", boundaryflag);
   AddDoubleProperty("launchPointX", launchPointX);
   AddDoubleProperty("launchPointY", launchPointY);
   AddDoubleProperty("launchPointZ", launchPointZ);
   AddDoubleProperty("focusPointX", focusPointX);
   AddDoubleProperty("focusPointY", focusPointY);
   AddDoubleProperty("focusPointZ", focusPointZ);
   AddDoubleProperty("trajectoryVectorX", trajectoryVectorX);
   AddDoubleProperty("trajectoryVectorY", trajectoryVectorY);
   AddDoubleProperty("trajectoryVectorZ", trajectoryVectorZ);
   AddDoubleProperty("radius", radius);
   AddDoubleProperty("waist", waist);
   AddDoubleProperty("sigma", m_Sigma);
   AddDoubleProperty("standardTissueAbsorption", initAbsorption);
   AddDoubleProperty("standardTissueScattering", initScattering);
   AddDoubleProperty("standardTissueAnisotropy", initAnisotropy);
   AddDoubleProperty("airThickness", m_ThicknessAir);
   AddDoubleProperty("skinThickness", m_ThicknessSkin);
 
   m_State = INITIALIZED;
 }
 
 void mitk::PhotoacousticVolume::AddDoubleProperty(std::string label, double value)
 {
   m_PropertyList->SetDoubleProperty(label.c_str(), value);
   mitk::CoreServices::GetPropertyPersistence()->AddInfo(label, mitk::PropertyPersistenceInfo::New(label));
 }
 
 void mitk::PhotoacousticVolume::AddIntProperty(std::string label, int value)
 {
   m_PropertyList->SetIntProperty(label.c_str(), value);
   mitk::CoreServices::GetPropertyPersistence()->AddInfo(label, mitk::PropertyPersistenceInfo::New(label));
 }
 
 mitk::Image::Pointer mitk::PhotoacousticVolume::ConvertToMitkImage()
 {
   AssertState(INITIALIZED);
 
   mitk::Image::Pointer resultImage = mitk::Image::New();
   mitk::PixelType TPixel = mitk::MakeScalarPixelType<double>();
   unsigned int* dimensionsOfImage = new unsigned int[4];
 
   // Copy dimensions
   dimensionsOfImage[0] = m_YDim;
   dimensionsOfImage[1] = m_XDim;
   dimensionsOfImage[2] = m_ZDim;
   dimensionsOfImage[3] = m_TDim;
 
   MITK_DEBUG << "Initializing image...";
   resultImage->Initialize(TPixel, 4, dimensionsOfImage, 1);
   MITK_DEBUG << "Initializing image...[Done]";
 
   mitk::Vector3D spacing;
   spacing.Fill(m_Spacing);
   resultImage->SetSpacing(spacing);
 
   MITK_DEBUG << "Set Import Volumes...";
   //Copy memory, deal with cleaning up memory yourself!
   resultImage->SetImportVolume(m_AbsorptionArray, 0, 0, mitk::Image::CopyMemory);
   resultImage->SetImportVolume(m_ScatteringArray, 1, 0, mitk::Image::CopyMemory);
   resultImage->SetImportVolume(m_AnisotropyArray, 2, 0, mitk::Image::CopyMemory);
   resultImage->SetImportVolume(m_SegmentationArray, 3, 0, mitk::Image::CopyMemory);
   MITK_DEBUG << "Set Import Volumes...[Done]";
 
   resultImage->SetPropertyList(m_PropertyList);
 
   return resultImage;
 }
 
 mitk::PhotoacousticVolume::PhotoacousticVolume()
 {
   m_State = UNINITIALIZED;
   m_XDim = 0;
   m_YDim = 0;
   m_ZDim = 0;
 
   m_AbsorptionAir = 0.0001;
   m_ScatteringAir = 1.0;
   m_AnisotropyAir = 1.0;
 
   m_Sigma = 0;
 
   m_PropertyList = mitk::PropertyList::New();
 }
 
 mitk::PhotoacousticVolume::~PhotoacousticVolume()
 {
   m_State = UNINITIALIZED;
   m_XDim = 0;
   m_YDim = 0;
   m_ZDim = 0;
 
   if(m_AbsorptionArray!=nullptr)
   {
     delete[] m_AbsorptionArray;
     m_AbsorptionArray = nullptr;
   }
   if(m_ScatteringArray!=nullptr)
   {
     delete[] m_ScatteringArray;
     m_ScatteringArray = nullptr;
   }
   if(m_AnisotropyArray!=nullptr)
   {
     delete[] m_AnisotropyArray;
     m_AnisotropyArray = nullptr;
   }
 
 }
 
 void mitk::PhotoacousticVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType)
 {
   AssertState(INITIALIZED);
   if(IsInsideVolume(x, y, z))
   {
     int index = z * m_XDim * m_YDim + x * m_YDim + y;
     m_AbsorptionArray[index] = absorption;
     m_ScatteringArray[index] = scattering;
     m_AnisotropyArray[index] = anisotropy;
     m_SegmentationArray[index] = segmentType;
   }
 }
 
 bool mitk::PhotoacousticVolume::IsInsideVolume(int x, int y, int z)
 {
   AssertState(INITIALIZED);
   return x >= 0 && x < m_XDim && y >= 0 && y < m_YDim && z >= 0 && z < m_ZDim;
 }
 
 double mitk::PhotoacousticVolume::GetVolumeAbsorptionValue(int x, int y, int z)
 {
   AssertState(INITIALIZED);
   return m_AbsorptionArray[z * m_XDim * m_YDim + x * m_YDim + y];
 }
 
 double mitk::PhotoacousticVolume::GetVolumeSegmentationValue(int x, int y, int z)
 {
   AssertState(INITIALIZED);
   return m_SegmentationArray[z * m_XDim * m_YDim + x * m_YDim + y];
 }
 
 void mitk::PhotoacousticVolume::FinalizeVolume(bool randomize, double percentage, bool useRngSeed, long rngSeed)
 {
   //Add air and skin areas
   double airvoxel = m_ThicknessAir / m_Spacing  / 10;
   double skinvoxel = 3 + m_ThicknessSkin / m_Spacing  / 10;
 
   for(int y = 0; y < m_YDim; y++)
   {
     for(int x = 0; x < m_XDim; x++)
     {
       for(int z = 0; z < airvoxel; z++)
       {
         int index = z * m_XDim * m_YDim + x * m_YDim + y;
         if(airvoxel-z < 0.9)
         {
           m_AbsorptionArray[index] = (1 - (airvoxel - z)) * m_AbsorptionArray[index] + (airvoxel - z) * m_AbsorptionAir;
           m_ScatteringArray[index] = (1 - (airvoxel - z)) * m_ScatteringArray[index] + (airvoxel - z) * m_ScatteringAir;
           m_AnisotropyArray[index] = (1 - (airvoxel - z)) * m_AnisotropyArray[index] + (airvoxel - z) * m_AnisotropyAir;
           if(airvoxel-z > 0.5)
             m_SegmentationArray[index] = SegmentationType::AIR;
         }
         else
         {
           m_AbsorptionArray[index] = m_AbsorptionAir;
           m_ScatteringArray[index] = m_ScatteringAir;
           m_AnisotropyArray[index] = m_AnisotropyAir;
           m_SegmentationArray[index] = SegmentationType::AIR;
         }
       }
 
       for(int z=airvoxel; z < skinvoxel; z++)
       {
         int index = z * m_XDim * m_YDim + x * m_YDim + y;
-        if(skinvoxel-z < 0.9)
+        if(skinvoxel-z < 1)
         {
           m_AbsorptionArray[index] = (1 - (skinvoxel - z)) * m_AbsorptionArray[index] + (skinvoxel - z) * m_AbsorptionSkin;
           m_ScatteringArray[index] = (1 - (skinvoxel - z)) * m_ScatteringArray[index] + (skinvoxel - z) * m_ScatteringSkin;
           m_AnisotropyArray[index] = (1 - (skinvoxel - z)) * m_AnisotropyArray[index] + (skinvoxel - z) * m_AnisotropySkin;
           if(skinvoxel-z > 0.5)
             m_SegmentationArray[index] = SegmentationType::SKIN;
         }
         else
         {
           m_AbsorptionArray[index] = m_AbsorptionSkin;
           m_ScatteringArray[index] = m_ScatteringSkin;
           m_AnisotropyArray[index] = m_AnisotropySkin;
           m_SegmentationArray[index] = SegmentationType::SKIN;
         }
       }
     }
   }
 
   //Randomize the coefficients if wanted
   if(!randomize)
     return;
 
   std::mt19937 rng;
   std::random_device randomDevice;
   if(useRngSeed)
   {
     rng.seed(rngSeed);
   }
   else
   {
     if(randomDevice.entropy()>0.1)
     {
       rng.seed(randomDevice());
     }
     else
     {
       rng.seed(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
     }
   }
   std::normal_distribution<> percentageDistribution(1, percentage/100);
 
   for(int y = 0; y < m_YDim; y++)
   {
     for(int x = 0; x < m_XDim; x++)
     {
       for(int z = 0; z < m_ZDim; z++)
       {
         int index = z * m_XDim * m_YDim + x * m_YDim + y;
         m_AbsorptionArray[index] = m_AbsorptionArray[index]*percentageDistribution(rng);
         m_ScatteringArray[index] = m_ScatteringArray[index]*percentageDistribution(rng);
       }
     }
   }
 
 }
 
 double mitk::PhotoacousticVolume::GetVolumeScatteringValue(int x, int y, int z)
 {
   AssertState(INITIALIZED);
   return m_ScatteringArray[z * m_XDim * m_YDim + x * m_YDim + y];
 }
 
 double mitk::PhotoacousticVolume::GetVolumeAnisotropyValue(int x, int y, int z)
 {
   AssertState(INITIALIZED);
   return m_AnisotropyArray[z * m_XDim * m_YDim + x * m_YDim + y];
 }
 
 void mitk::PhotoacousticVolume::GaussianBlur3D()
 {
   GaussianBlur3D(m_AbsorptionArray, m_YDim, m_XDim, m_ZDim, m_Sigma);
   GaussianBlur3D(m_ScatteringArray, m_YDim, m_XDim, m_ZDim, m_Sigma);
   GaussianBlur3D(m_AnisotropyArray, m_YDim, m_XDim, m_ZDim, m_Sigma);
 }
 
 /**
  * @brief Fast 3D Gaussian convolution IIR approximation
  * @param volume
  * @param width
  * @param height
  * @param depth
  * @param sigma
  * @author Pascal Getreuer <getreuer@gmail.com>
  *
  * Copyright (c) 2011, Pascal Getreuer
  * All rights reserved.
  *
  * This program is free software: you can redistribute it and/or modify it
  * under the terms of the simplified BSD license.
  *
  * You should have received a copy of these licenses along with this program.
  * If not, see <http://www.opensource.org/licenses/bsd-license.html>.
  */
 void mitk::PhotoacousticVolume::GaussianBlur3D(double *volume, long width, long height, long depth,
                                                double sigma)
 {
   const long plane = width*height;
   const long numel = plane*depth;
   double lambda, dnu;
   double nu, boundaryscale, postscale;
   double *ptr;
   long i, x, y, z;
   int step;
 
   if(sigma <= 0)
     return;
 
   lambda = (sigma*sigma)/(8.0);
   dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda))/(2.0*lambda);
   nu = dnu;
   boundaryscale = 1.0/(1.0 - dnu);
   postscale = pow(dnu/lambda,12);
 
   /* Filter horizontally along each row */
   for(z = 0; z < depth; z++)
   {
     for(y = 0; y < height; y++)
     {
       for(step = 0; step < 4; step++)
       {
         ptr = volume + width*(y + height*z);
         ptr[0] *= boundaryscale;
 
         /* Filter rightwards */
         for(x = 1; x < width; x++)
         {
           ptr[x] += nu*ptr[x - 1];
         }
 
         ptr[x = width - 1] *= boundaryscale;
             /* Filter leftwards */
             for(; x > 0; x--)
         {
           ptr[x - 1] += nu*ptr[x];
         }
       }
     }
   }
   /* Filter vertically along each column */
   for(z = 0; z < depth; z++)
   {
     for(x = 0; x < width; x++)
     {
       for(step = 0; step < 4; step++)
       {
         ptr = volume + x + plane*z;
         ptr[0] *= boundaryscale;
 
         /* Filter downwards */
         for(i = width; i < plane; i += width)
         {
           ptr[i] += nu*ptr[i - width];
         }
 
         ptr[i = plane - width] *= boundaryscale;
 
             /* Filter upwards */
             for(; i > 0; i -= width)
         {
           ptr[i - width] += nu*ptr[i];
         }
       }
     }
   }
 
   /* Filter along z-dimension */
   for(y = 0; y < height; y++)
   {
     for(x = 0; x < width; x++)
     {
       for(step = 0; step < 4; step++)
       {
         ptr = volume + x + width*y;
         ptr[0] *= boundaryscale;
 
         for(i = plane; i < numel; i += plane)
         {
           ptr[i] += nu*ptr[i - plane];
         }
 
         ptr[i = numel - plane] *= boundaryscale;
 
             for(; i > 0; i -= plane)
         {
           ptr[i - plane] += nu*ptr[i];
         }
       }
     }
   }
 
   for(i = 0; i < numel; i++)
   {
     volume[i] *= postscale;
   }
 
   return;
 }
diff --git a/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp
index d2f272831c..c33a9f6a0c 100644
--- a/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp
+++ b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp
@@ -1,140 +1,139 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #include "mitkPhotoacousticSmartVector.h"
 #include "chrono"
-#include "itkAffineTransform.h"
-
+#include "math.h"
 mitk::PhotoacousticSmartVector::PhotoacousticSmartVector()
 {
   srand(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
   m_Vector.Fill(0);
 }
 
 mitk::PhotoacousticSmartVector::~PhotoacousticSmartVector()
 {
   m_Vector.Fill(0);
 }
 
 double mitk::PhotoacousticSmartVector::GetNorm()
 {
   return m_Vector.GetNorm();
 }
 
 double mitk::PhotoacousticSmartVector::GetElement(unsigned short index)
 {
   return m_Vector.GetElement(index);
 }
 
 void mitk::PhotoacousticSmartVector::SetElement(unsigned short index, double value)
 {
   m_Vector.SetElement(index, value);
 }
 
 void mitk::PhotoacousticSmartVector::Normalize()
 {
   double norm = m_Vector.GetNorm();
   m_Vector.SetElement(0, m_Vector.GetElement(0)/norm);
   m_Vector.SetElement(1, m_Vector.GetElement(1)/norm);
   m_Vector.SetElement(2, m_Vector.GetElement(2)/norm);
 }
 
 void mitk::PhotoacousticSmartVector::SetValue(mitk::PhotoacousticSmartVector::Pointer value)
 {
   m_Vector.SetElement(0, value->GetElement(0));
   m_Vector.SetElement(1, value->GetElement(1));
   m_Vector.SetElement(2, value->GetElement(2));
 }
 
-void mitk::PhotoacousticSmartVector::RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937 rng)
+void mitk::PhotoacousticSmartVector::RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937* rng)
 {
   std::uniform_real_distribution<> range(-percentage, percentage);
-  m_Vector.SetElement(0, m_Vector.GetElement(0) + (bendingFactor * range(rng)));
-  m_Vector.SetElement(1, m_Vector.GetElement(1) + (bendingFactor * range(rng)));
-  m_Vector.SetElement(2, m_Vector.GetElement(2) + (bendingFactor * range(rng)));
+  m_Vector.SetElement(0, m_Vector.GetElement(0) + (bendingFactor * range(*rng)));
+  m_Vector.SetElement(1, m_Vector.GetElement(1) + (bendingFactor * range(*rng)));
+  m_Vector.SetElement(2, m_Vector.GetElement(2) + (bendingFactor * range(*rng)));
 }
 
-void mitk::PhotoacousticSmartVector::Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937 rng)
+void mitk::PhotoacousticSmartVector::Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937* rng)
 {
   std::uniform_real_distribution<> rangeX(xLowerLimit, xUpperLimit);
   std::uniform_real_distribution<> rangeY(yLowerLimit, yUpperLimit);
   std::uniform_real_distribution<> rangeZ(zLowerLimit, zUpperLimit);
-  m_Vector.SetElement(0, rangeX(rng));
-  m_Vector.SetElement(1, rangeY(rng));
-  m_Vector.SetElement(2, rangeZ(rng));
+  m_Vector.SetElement(0, rangeX(*rng));
+  m_Vector.SetElement(1, rangeY(*rng));
+  m_Vector.SetElement(2, rangeZ(*rng));
 }
 
-void mitk::PhotoacousticSmartVector::Randomize(double xLimit, double yLimit, double zLimit, std::mt19937 rng)
+void mitk::PhotoacousticSmartVector::Randomize(double xLimit, double yLimit, double zLimit, std::mt19937* rng)
 {
   Randomize(0, xLimit, 0, yLimit, 0, zLimit, rng);
 }
 
-void mitk::PhotoacousticSmartVector::Randomize(std::mt19937 rng)
+void mitk::PhotoacousticSmartVector::Randomize(std::mt19937* rng)
 {
   Randomize(-1, 1, -1, 1, -1, 1, rng);
 }
 
 void mitk::PhotoacousticSmartVector::PrintSelf(std::ostream& os, itk::Indent indent) const
 {
   os << "X: " << m_Vector.GetElement(0) << std::endl;
   os << "Y: " << m_Vector.GetElement(1) << std::endl;
   os << "Z: " << m_Vector.GetElement(2) << std::endl;
   os << "Length: " << m_Vector.GetNorm() << std::endl;
 }
 
-void mitk::PhotoacousticSmartVector::Rotate(double xAngle, double yAngle, double zAngle)
+void mitk::PhotoacousticSmartVector::Rotate(double thetaChange, double phiChange)
 {
-  itk::AffineTransform<double, 3>::Pointer transform = itk::AffineTransform<double,3>::New();
+  MITK_INFO << "Vector before rotation: ("<<GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")";
+  if(thetaChange==0 && phiChange == 0)
+    return;
 
-  itk::Vector<double, 3> xAxis;
-  xAxis.Fill(0);
-  xAxis[0] = 1;
+  double x = GetElement(0);
+  double y = GetElement(1);
+  double z = GetElement(2);
 
-  itk::Vector<double, 3> yAxis;
-  yAxis.Fill(0);
-  yAxis[1] = 1;
+  double r = sqrt(x*x + y*y + z*z);
+  if(r==0)
+    return;
 
-  itk::Vector<double, 3> zAxis;
-  zAxis.Fill(0);
-  zAxis[2] = 1;
+  double theta = acos(z/r);
+  double phi = atan2(y, x);
 
-  transform->Rotate3D(xAxis, xAngle);
-  transform->Rotate3D(yAxis, yAngle);
-  transform->Rotate3D(zAxis, zAngle);
+  theta += thetaChange;
+  phi += phiChange;
 
-  itk::Vector<double,3> rotatedVector = transform->TransformVector(m_Vector);
+  SetElement(0, r * sin(theta) * cos(phi));
+  SetElement(1, r * sin(theta) * sin(phi));
+  SetElement(2, r * cos(theta));
 
-  m_Vector.SetElement(0, rotatedVector.GetElement(0));
-  m_Vector.SetElement(1, rotatedVector.GetElement(1));
-  m_Vector.SetElement(2, rotatedVector.GetElement(2));
+  MITK_INFO << "Vector after rotation: ("<<GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")";
 
 }
 
 void mitk::PhotoacousticSmartVector::Scale(double factor)
 {
   m_Vector.SetElement(0, m_Vector.GetElement(0)*factor);
   m_Vector.SetElement(1, m_Vector.GetElement(1)*factor);
   m_Vector.SetElement(2, m_Vector.GetElement(2)*factor);
 }
 
 mitk::PhotoacousticSmartVector::Pointer mitk::PhotoacousticSmartVector::Clone()
 {
   mitk::PhotoacousticSmartVector::Pointer returnVector = mitk::PhotoacousticSmartVector::New();
   returnVector->SetElement(0, this->GetElement(0));
   returnVector->SetElement(1, this->GetElement(1));
   returnVector->SetElement(2, this->GetElement(2));
   return returnVector;
 }
diff --git a/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h
index 18f0df9550..87b16d2b8b 100644
--- a/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h
+++ b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h
@@ -1,125 +1,124 @@
 /*===================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center,
 Division of Medical and Biological Informatics.
 All rights reserved.
 
 This software is distributed WITHOUT ANY WARRANTY; without
 even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A PARTICULAR PURPOSE.
 
 See LICENSE.txt or http://www.mitk.org for details.
 
 ===================================================================*/
 
 #ifndef MITKSMARTVECTOR_H
 #define MITKSMARTVECTOR_H
 
 #include <mitkVector.h>
 #include <random>
 
 #include <MitkPhotoacousticSimulationExports.h>
 
 //Includes for smart pointer usage
 #include "mitkCommon.h"
 #include "itkLightObject.h"
 
 namespace mitk {
 
 
 class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticSmartVector : public itk::LightObject
 {
 public:
   mitkClassMacroItkParent(mitk::PhotoacousticSmartVector, itk::LightObject)
   itkFactorylessNewMacro(Self)
 
   /**
    * @brief GetNorm calculates the length of this vector.
    * @return the euclidean norm
    */
   double GetNorm();
 
   double GetElement(unsigned short index);
   void SetElement(unsigned short index, double value);
 
   /**
    * @brief Normalize normalizes this vector. After calling this GetNorm() will return 1.
    */
   void Normalize();
   void SetValue(mitk::PhotoacousticSmartVector::Pointer value);
 
   /**
    * @brief RandomizeByPercentage alters this vector randomly by [-percentage, percentage] of the bendingFactor.
    *
    * @param percentage
    * @param bendingFactor
    */
-  void RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937 rng);
+  void RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937* rng);
 
   /**
    * @brief Randomize randomizes this vector to be [lowerLimit, upperLimit] in each element
    *
    * @param xLowerLimit
    * @param xUpperLimit
    * @param yLowerLimit
    * @param yUpperLimit
    * @param zLowerLimit
    * @param zUpperLimit
    */
-  void Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937 rng);
+  void Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937* rng);
 
   /**
    * @brief Randomize randomizes this vector to be [0, limit] in each element
    *
    * @param xLimit
    * @param yLimit
    * @param zLimit
    */
-  void Randomize(double xLimit, double yLimit, double zLimit, std::mt19937 rng);
+  void Randomize(double xLimit, double yLimit, double zLimit, std::mt19937* rng);
 
   /**
    * @brief Randomize randomizes this vector to be [-1, 1] in each element
    */
-  void Randomize(std::mt19937 rng);
+  void Randomize(std::mt19937* rng);
 
   /**
    * @brief Rotate rotates this PhotoacousticSmartVector around the x, y and z axis with the given angles in radians
    *
-   * @param xAngle x angle in radians
-   * @param yAngle y angle in radians
-   * @param zAngle z angle in radians
+   * @param thetaChange rotation of the inclination angle in radians
+   * @param phiChange rotation of the azimuthal angle in radians
    */
-  void Rotate(double xAngle, double yAngle, double zAngle);
+  void Rotate(double xAngle, double yAngle);
 
   /**
    * @brief Scale scales this PhotoacousticSmartVector with the given factor
    *
    * @param factor the scaling factor
    *
    * If a negative number is provided, the direction of the vector will be inverted.
    */
   void Scale(double factor);
 
   /**
    * @brief Clone create a deep copy of this vector
    *
    * @return a new vector with the same values.
    */
   mitk::PhotoacousticSmartVector::Pointer Clone();
 
 
 protected:
   PhotoacousticSmartVector();
   virtual ~PhotoacousticSmartVector();
 
   virtual void PrintSelf(std::ostream& os, itk::Indent indent) const override;
 
 private:
   mitk::Vector3D m_Vector;
 };
 
 }
 
 #endif // MITKSMARTVECTOR_H