diff --git a/code/algorithms/files.cmake b/code/algorithms/files.cmake
index 458f698..981351b 100644
--- a/code/algorithms/files.cmake
+++ b/code/algorithms/files.cmake
@@ -1,12 +1,17 @@
 SET(CPP_FILES 
   rttbDoseStatistics.cpp
+  rttbDoseStatisticsCalculator.cpp
   rttbArithmetic.cpp
   )
 
 SET(H_FILES 
   rttbDoseStatistics.h
+  rttbDoseStatisticsCalculator.h
   rttbArithmetic.h
-  rttbArithmetic.tpp
   rttbBinaryFunctorDoseAccessor.h
-  rttbBinaryFunctorDoseAccessor.tpp
+)
+
+SET(TXX_FILES
+	rttbArithmetic.tpp
+	rttbBinaryFunctorDoseAccessor.tpp
 )
diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp
index ac63f89..14fe157 100644
--- a/code/algorithms/rttbDoseStatistics.cpp
+++ b/code/algorithms/rttbDoseStatistics.cpp
@@ -1,453 +1,314 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE.  See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date    $Date$ (last change date)
 // @author  $Author$ (last changed by)
 */
 
-#include <string>
-#include <vector>
-#include <exception>
-#include <math.h>
-#include <deque>
-#include <algorithm>
-#include <iostream>
-#include <map>
-
 #include "rttbDoseStatistics.h"
-#include "rttbNullPointerException.h"
-#include "rttbInvalidDoseException.h"
+#include "rttbDataNotAvailableException.h"
 
 namespace rttb
 {
 
 	namespace algorithms
 	{
-		DoseStatistics::DoseStatistics()
+		DoseStatistics::DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean,
+		                               DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume,
+		                               ResultListPointer maximumVoxelPositions /*= ResultListPointer()*/,
+		                               ResultListPointer minimumVoxelPositions /*= ResultListPointer()*/,
+		                               DoseToVolumeFunctionType Dx /*= std::map<DoseTypeGy, VolumeType>()*/,
+		                               VolumeToDoseFunctionType Vx /*= std::map<VolumeType, DoseTypeGy>()*/,
+		                               VolumeToDoseFunctionType MOHx /*= std::map<VolumeType, DoseTypeGy>()*/,
+		                               VolumeToDoseFunctionType MOCx /*= std::map<VolumeType, DoseTypeGy>()*/,
+		                               VolumeToDoseFunctionType MaxOHx /*= std::map<VolumeType, DoseTypeGy>()*/,
+		                               VolumeToDoseFunctionType MinOCx /*= std::map<VolumeType, DoseTypeGy>()*/) : _minimum(minimum),
+			_maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume),
+			_maximumVoxelPositions(maximumVoxelPositions), _minimumVoxelPositions(minimumVoxelPositions), _Dx(Dx), _Vx(Vx),
+			_MOHx(MOHx),
+			_MOCx(MOCx), _MaxOHx(MaxOHx), _MinOCx(MinOCx)
 		{
-			initSuccess = false;
+
 		}
 
-		DoseStatistics::DoseStatistics(DoseIteratorPointer aDoseIterator)
+
+		DoseStatistics::~DoseStatistics()
 		{
-			_doseIterator = aDoseIterator;
-			initSuccess = false;
-			this->init();
+
 		}
 
-		void DoseStatistics::setDoseIterator(DoseIteratorPointer aDoseIterator)
+
+		void DoseStatistics::setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions)
 		{
-			_doseIterator = aDoseIterator;
-			initSuccess = false;
-			this->init();
+			_minimumVoxelPositions = minimumVoxelPositions;
 		}
 
-		DoseStatistics::DoseIteratorPointer DoseStatistics::getDoseIterator() const
+		void DoseStatistics::setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions)
 		{
-			return _doseIterator;
+			_maximumVoxelPositions = maximumVoxelPositions;
 		}
 
-
-		bool DoseStatistics::init()
+		void DoseStatistics::setDx(const DoseToVolumeFunctionType& DxValues)
 		{
-
-			if (!_doseIterator)
-			{
-				throw core::NullPointerException("_doseIterator must not be NULL!");
-			}
-
-			doseVector.clear();
-			voxelProportionVector.clear();
-
-			std::multimap<double, int> doseValueVSIndexMap;
-			std::vector<double> voxelProportionVectorTemp;
-
-
-			this->_maximum = 0;
-			this->_mean = 0;
-			this->_stdDeviation = 0;
-			this->_variance = 0;
-
-			float sum = 0;
-			_numVoxels = 0;
-			float squareSum = 0;
-
-			_doseIterator->reset();
-			int i = 0;
-			DoseTypeGy doseValue = 0;
-
-			while (_doseIterator->isPositionValid())
-			{
-				doseValue = _doseIterator->getCurrentDoseValue();
-
-				if (i == 0)
-				{
-					_minimum = doseValue;
-				}
-
-				rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction();
-				sum += doseValue * voxelProportion;
-				_numVoxels += voxelProportion;
-				squareSum += doseValue * doseValue * voxelProportion;
-
-				if (doseValue > this->_maximum)
-				{
-					_maximum = doseValue;
-				}
-				else if (doseValue < this->_minimum)
-				{
-					_minimum = doseValue;
-				}
-
-				voxelProportionVectorTemp.push_back(voxelProportion);
-				doseValueVSIndexMap.insert(std::pair<double, int>(doseValue, i));
-
-				i++;
-				_doseIterator->next();
-			}
-
-			if (_numVoxels != 0)
-			{
-				_mean = sum / _numVoxels;
-				_variance = (squareSum / _numVoxels - _mean * _mean);
-
-				if (_variance < errorConstant)
-				{
-					_stdDeviation = 0;
-				}
-				else
-				{
-					_stdDeviation = pow(_variance, 0.5);
-				}
-			}
-
-			//sort dose values and corresponding volume fractions in member variables
-			std::multimap<double, int>::iterator it;
-
-			for (it = doseValueVSIndexMap.begin(); it != doseValueVSIndexMap.end(); ++it)
-			{
-				doseVector.push_back((float)(*it).first);
-				voxelProportionVector.push_back(voxelProportionVectorTemp.at((*it).second));
-			}
-
-			initSuccess = true;
-
-			return true;
+			_Dx = DxValues;
 		}
 
-		double DoseStatistics::getNumberOfVoxels()
+		void DoseStatistics::setVx(const VolumeToDoseFunctionType& VxValues)
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
-			}
+			_Vx = VxValues;
+		}
 
-			return _numVoxels;
+		void DoseStatistics::setMOHx(const VolumeToDoseFunctionType& MOHxValues)
+		{
+			_MOHx = MOHxValues;
 		}
 
-		DoseStatisticType DoseStatistics::getMaximum(ResultListPointer maxVoxelVector) const
+		void DoseStatistics::setMOCx(const VolumeToDoseFunctionType& MOCxValues)
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
+			_MOCx = MOCxValues;
+		}
 
-			}
+		void DoseStatistics::setMaxOHx(const VolumeToDoseFunctionType& MaxOHValues)
+		{
+			_MaxOHx = MaxOHValues;
+		}
 
-			if (maxVoxelVector == NULL)
-			{
-				throw core::NullPointerException("resultsVector must not be NULL! ");
-			}
+		void DoseStatistics::setMinOCx(const VolumeToDoseFunctionType& MinOCValues)
+		{
+			_MinOCx = MinOCValues;
+		}
 
-			if (maxVoxelVector->size() == 0)
-			{
-				this->_doseIterator->reset();
-				DoseTypeGy doseValue = 0;
+		unsigned int DoseStatistics::getNumberOfVoxels() const
+		{
+			return _numVoxels;
+		}
 
-				while (_doseIterator->isPositionValid())
-				{
-					doseValue = _doseIterator->getCurrentDoseValue();
 
-					if (doseValue == _maximum)
-					{
-						VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID();
-						std::pair<DoseTypeGy, VoxelGridID> voxel(doseValue, currentID);
-						maxVoxelVector->push_back(voxel);
-					}
+		VolumeType DoseStatistics::getVolume() const
+		{
+			return _volume;
+		}
 
-					_doseIterator->next();
-				}
-			}
 
+		DoseStatisticType DoseStatistics::getMaximum() const
+		{
 			return _maximum;
 		}
 
-		DoseStatisticType DoseStatistics::getMinimum(ResultListPointer minVoxelVector, int number) const
+		DoseStatisticType DoseStatistics::getMinimum() const
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
-
-			}
-
-			if (minVoxelVector == NULL)
-			{
-				throw core::NullPointerException("resultsVector must not be NULL! ");
-			}
-
-			/*! @todo: Architecture Annotation:
-				Finding the positions for the minimum only once reduces computation time,
-				but will require sensible use by the programmers. To be save the output vector
-				minVoxelVector will be always cleared here to garantee that no false values are
-				presented. This change may be revoced to increase computation speed later on
-				(only compute if(minVoxelVector->size()==0)).
-			*/
-			minVoxelVector->clear();
-			int count = 0;
-			this->_doseIterator->reset();
-			DoseTypeGy doseValue = 0;
-
-			while (_doseIterator->isPositionValid() && count < number)
-			{
-				doseValue = _doseIterator->getCurrentDoseValue();
-
-				if (doseValue == _minimum)
-				{
-					VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID();
-					std::pair<DoseTypeGy, VoxelGridID> voxel(doseValue, currentID);
-					minVoxelVector->push_back(voxel);
-					count++;
-				}
-
-				_doseIterator->next();
-			}
-
 			return _minimum;
 		}
 
 		DoseStatisticType DoseStatistics::getMean() const
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
-			}
-
 			return _mean;
 		}
 
 		DoseStatisticType DoseStatistics::getStdDeviation() const
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
-			}
-
 			return _stdDeviation;
 		}
 
 		DoseStatisticType DoseStatistics::getVariance() const
 		{
-			if (!initSuccess)
-			{
-				throw core::InvalidDoseException("DoseStatistics is not initialized: set dose using setDoseIterator()! ");
-			}
-
-			return _variance;
+			return _stdDeviation * _stdDeviation;
 		}
 
-		VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute) const
+		VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue,
+		                                 DoseTypeGy& nearestXDose) const
 		{
-			rttb::FractionType count = 0;
-			_doseIterator->reset();
-
-			DoseTypeGy currentDose = 0;
+			return getValue(_Vx, xDoseAbsolute, findNearestValue, nearestXDose);
 
-			while (_doseIterator->isPositionValid())
-			{
-				currentDose = _doseIterator->getCurrentDoseValue();
+		}
 
-				if (currentDose >= xDoseAbsolute)
-				{
-					count += _doseIterator->getCurrentRelevantVolumeFraction();
-				}
+		VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute) const
+		{
+			DoseTypeGy dummy;
+			return getValue(_Vx, xDoseAbsolute, false, dummy);
+		}
 
-				_doseIterator->next();
-			}
+		DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute, bool findNearestValue,
+		                                 VolumeType& nearestXVolume) const
+		{
+			return getValue(_Dx, xVolumeAbsolute, findNearestValue, nearestXVolume);
+		}
 
-			return count * this->_doseIterator->getCurrentVoxelVolume();
+		DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute) const
+		{
+			VolumeType dummy;
+			return getValue(_Dx, xVolumeAbsolute, false, dummy);
 		}
 
-		DoseTypeGy DoseStatistics::getDx(DoseTypeGy xVolumeAbsolute) const
+		DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue,
+		                                   VolumeType& nearestXVolume) const
 		{
-			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
-			DoseTypeGy resultDose = 0;
+			return getValue(_MOHx, xVolumeAbsolute, findNearestValue, nearestXVolume);
+		}
 
-			double countVoxels = 0;
-			int i = doseVector.size() - 1;
+		DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute) const
+		{
+			VolumeType dummy;
+			return getValue(_MOHx, xVolumeAbsolute, false, dummy);
+		}
 
-			for (; i >= 0; i--)
-			{
-				countVoxels += voxelProportionVector.at(i);
+		DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue,
+		                                   VolumeType& maximumDistanceFromOriginalVolume) const
+		{
+			return getValue(_MOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume);
+		}
 
-				if (countVoxels >= noOfVoxel)
-				{
-					break;
-				}
-			}
+		DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute) const
+		{
+			VolumeType dummy;
+			return getValue(_MOCx, xVolumeAbsolute, false, dummy);
+		}
 
-			if (i >= 0)
-			{
-				resultDose = doseVector.at(i);
-			}
-			else
-			{
-				resultDose = _minimum;
-			}
+		DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue,
+		                                     VolumeType& maximumDistanceFromOriginalVolume) const
+		{
+			return getValue(_MaxOHx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume);
+		}
 
-			return resultDose;
+		DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute) const
+		{
+			VolumeType dummy;
+			return getValue(_MaxOHx, xVolumeAbsolute, false, dummy);
 		}
 
-		DoseTypeGy DoseStatistics::getMOHx(DoseTypeGy xVolumeAbsolute) const
+		DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue,
+		                                     VolumeType& maximumDistanceFromOriginalVolume) const
 		{
-			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+			return getValue(_MinOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume);
+		}
 
+		DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute) const
+		{
+			VolumeType dummy;
+			return getValue(_MinOCx, xVolumeAbsolute, false, dummy);
+		}
 
-			if (noOfVoxel == 0)
+		double DoseStatistics::getValue(std::map<double, double> aMap, double key,
+		                                bool findNearestValueInstead, double& storedKey) const
+		{
+			if (aMap.find(key) != aMap.end())
 			{
-				return 0;
+				return aMap.find(key)->second;
 			}
 			else
 			{
-				double countVoxels = 0;
-				double sum = 0;
-
-				for (int i = doseVector.size() - 1; i >= 0; i--)
+				//value not in map. We have to find the nearest value
+				if (aMap.empty())
 				{
-					double voxelProportion = voxelProportionVector.at(i);
-					countVoxels += voxelProportion;
-					sum += doseVector.at(i) * voxelProportion;
-
-					if (countVoxels >= noOfVoxel)
+					throw core::DataNotAvailableException("No Vx values are defined");
+				}
+				else
+				{
+					if (findNearestValueInstead)
 					{
-						break;
+						auto iterator = findNearestKeyInMap(aMap, key);
+						storedKey = iterator->first;
+						return iterator->second;
+					}
+					else
+					{
+						throw core::DataNotAvailableException("No Vx value with required dose is defined");
 					}
 				}
-
-				return (DoseTypeGy)(sum / noOfVoxel);
 			}
 		}
 
-		DoseTypeGy DoseStatistics::getMOCx(DoseTypeGy xVolumeAbsolute) const
+		std::map<double, double>::const_iterator DoseStatistics::findNearestKeyInMap(std::map<double, double> aMap,
+		        double key) const
 		{
-			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+			double minDistance = 1e19;
+			double minDistanceLast = 1e20;
 
+			auto iterator = std::begin(aMap);
 
-			if (noOfVoxel == 0)
-			{
-				return 0;
-			}
-			else
+			while (iterator != std::end(aMap))
 			{
-				double countVoxels = 0;
-				double sum = 0;
-				std::vector<DoseTypeGy>::const_iterator it = doseVector.begin();
-				std::vector<double>::const_iterator itD = voxelProportionVector.begin();
+				minDistanceLast = minDistance;
+				minDistance = fabs(iterator->first - key);
 
-				for (; it != doseVector.end(); ++it, ++itD)
+				if (minDistanceLast > minDistance)
 				{
-					double voxelProportion = *itD;
-					countVoxels += voxelProportion;
-					sum += (*it) * voxelProportion;
-
-					if (countVoxels >= noOfVoxel)
+					++iterator;
+				}
+				else
+				{
+					if (iterator != std::begin(aMap))
+					{
+						--iterator;
+						return iterator;
+					}
+					else
 					{
-						break;
+						return std::begin(aMap);
 					}
 				}
-
-				return (DoseTypeGy)(sum / noOfVoxel);
 			}
+
+			--iterator;
+			return iterator;
 		}
 
-		DoseTypeGy DoseStatistics::getMaxOHx(DoseTypeGy xVolumeAbsolute) const
+		DoseStatistics::ResultListPointer DoseStatistics::getMaximumPositions() const
 		{
-			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
-			DoseTypeGy resultDose = 0;
-
-			double countVoxels = 0;
-			int i = doseVector.size() - 1;
-
-			for (; i >= 0; i--)
-			{
-				countVoxels += voxelProportionVector.at(i);
-
-				if (countVoxels >= noOfVoxel)
-				{
-					break;
-				}
-			}
+			return _maximumVoxelPositions;
+		}
 
-			if (i - 1 >= 0)
-			{
-				resultDose = doseVector.at(i - 1);
-			}
+		DoseStatistics::ResultListPointer DoseStatistics::getMinimumPositions() const
+		{
+			return _minimumVoxelPositions;
+		}
 
-			return resultDose;
+		DoseStatistics::DoseToVolumeFunctionType DoseStatistics::getAllVx() const
+		{
+			return _Vx;
 		}
 
-		DoseTypeGy DoseStatistics::getMinOCx(DoseTypeGy xVolumeAbsolute) const
+		DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllDx() const
 		{
-			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
-			DoseTypeGy resultDose = 0;
+			return _Dx;
+		}
 
-			double countVoxels = 0;
-			std::vector<DoseTypeGy>::const_iterator it = doseVector.begin();
-			std::vector<double>::const_iterator itD = voxelProportionVector.begin();
+		DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMOHx() const
+		{
+			return _MOHx;
+		}
 
-			for (; itD != voxelProportionVector.end(); ++itD, ++it)
-			{
-				countVoxels += *itD;
+		DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMOCx() const
+		{
+			return _MOCx;
+		}
 
-				if (countVoxels >= noOfVoxel)
-				{
-					break;
-				}
-			}
+		DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMaxOHx() const
+		{
+			return _MaxOHx;
+		}
 
-			if (it != doseVector.end())
-			{
-				++it;
+		DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMinOCx() const
+		{
+			return _MinOCx;
+		}
 
-				if (it != doseVector.end())
-				{
-					resultDose = *it;
-				}
-				else
-				{
-					resultDose = (DoseTypeGy)_maximum;
-				}
-			}
-			else
-			{
-				resultDose = (DoseTypeGy)_maximum;
-			}
 
-			return resultDose;
-		}
 
 	}//end namespace algorithms
 }//end namespace rttb
 
diff --git a/code/algorithms/rttbDoseStatistics.h b/code/algorithms/rttbDoseStatistics.h
index cc7de52..5a497cc 100644
--- a/code/algorithms/rttbDoseStatistics.h
+++ b/code/algorithms/rttbDoseStatistics.h
@@ -1,164 +1,205 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE.  See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date    $Date$ (last change date)
 // @author  $Author$ (last changed by)
 */
 
 #ifndef __DOSE_STATISTICS_H
 #define __DOSE_STATISTICS_H
 
+#include <boost/make_shared.hpp>
 #include <vector>
+#include <map>
 
 #include "rttbBaseType.h"
 #include "rttbDoseIteratorInterface.h"
 
-namespace rttb{
+namespace rttb
+{
 
-	namespace algorithms{
+	namespace algorithms
+	{
 
 		/*! @class DoseStatistics
-		@brief This is a class calculating different statistical values from a rt dose distribution
-		These values range fram standard statistical values such as minimum, maximum and mean to more 
-		dose specific properties such as Vx (volume iradiated with a dose >=x), Dx (minimal dose delivered 
-		to x% of the VOI, and MOHx (mean in the hottest volume).
+		@brief This is a data class storing different statistical values from a rt dose distribution
+		@sa DoseStatisticsCalculator
 		*/
 		class DoseStatistics
-			{
-			public:
-				typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer;
-				typedef boost::shared_ptr<std::vector<std::pair<DoseTypeGy,VoxelGridID> > > ResultListPointer;
-
-			private: 
-				DoseIteratorPointer _doseIterator;
-				DoseStatisticType _maximum;
-				DoseStatisticType _minimum;
-				DoseStatisticType _mean;
-				DoseStatisticType _stdDeviation;
-				DoseStatisticType _variance;
-				DoseStatisticType _numVoxels;
-
-				/* Contains relevant dose values sorted in descending order.
-				*/
-				std::vector<DoseTypeGy> doseVector;
-
-				/*! Contains the corresponding voxel proprtions to the values in doseVector. 
-				*/
-				std::vector<double> voxelProportionVector;
-
-				bool initSuccess;
-
-				/*! @brief Calculation of basic dose statistics. It will be called in Constructor.
-					@exception NullPointerException Thrown if _doseIterator is NULL.
-				*/
-				bool init();
-
-
-
-			public:
-				/*! @brief Standard Constructor
-				*/
-				DoseStatistics();
-
-				/*! @brief Constructor
-				*/
-				DoseStatistics(DoseIteratorPointer aDoseIterator);
-
-				/*! @brief Set new dose data for statistics. Statistics will be re-initialized.
-				*/
-				void setDoseIterator(DoseIteratorPointer aDoseIterator);
-
-				/*! @brief Get dose iterator.
-				*/
-				DoseIteratorPointer getDoseIterator() const;
-
-
-				/*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy.
-        @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getNumberOfVoxels();
-
-				/*! @brief Get the maximum of the current dose distribution.
-					@param maxVoxelVector: vector of all voxel with the dose=maximum.
-					@return Return the maximum dose in Gy, or -1 if an error occured.
-          @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getMaximum( ResultListPointer maxVoxelVector) const;
-
-				/*! @brief Get the minimum of the current dose distribution.
-					@param minVoxelVector: vector of all voxel with the dose=minimum.
-					@param number: the number of dose voxels with dose=minimum that are stored in minVoxelVector.
-					Only the first occurences are stored. The default value is 100.
-					@return Return the minimum dose in Gy, or -1 if an error occured.
-          @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getMinimum( ResultListPointer minVoxelVector, int number=100) const;
-
-				/*! @brief Get the mean of the current dose distribution.
-					@return Return the mean dose in Gy, or -1 if an error occured.
-          @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getMean() const;
-
-				/*! @brief Get the standard deviation of the current dose distribution.
-					@return Return the standard deviation in Gy, or -1 if an error occured.
-          @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getStdDeviation() const;
-
-				/*! @brief Get the variance of of the current dose distribution.
-					@return Return the variance in Gy, or -1 if an error occured.
-          @exception InvalidDoseException Thrown if statistic was not initialized.
-				*/
-				DoseStatisticType getVariance() const;
-
-				/*! @brief Get Vx: the volume irradiated with a dose >= x.
-					@return Return absolute volume in absolute cm3.
-				*/
-				VolumeType getVx(DoseTypeGy xDoseAbsolute) const;
-
-				/*! @brief Get Dx: the minimal dose delivered to part x of the current volume.
-					@return Return dose value in Gy.
-				*/	
-				DoseTypeGy getDx(VolumeType xVolumeAbsolute) const;
-
-				/*! @brief Get MOHx: mean dose of the hottest x voxels.
-					@return Return dose value in Gy.
-				*/	
-				DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) const;
-
-				/*! @brief Get MOCx: mean dose of the coldest x voxels.
-					@return Return dose value in Gy.
-				*/
-				DoseTypeGy getMOCx(VolumeType xVolumeAbsolute) const;
-
-				/*! @brief Get MaxOHx: Maximum outside of the hottest x voxels.
-					@return Return dose value in Gy.
-				*/
-				DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute) const;
-
-				/*! @brief Get MinOCx: Minimum outside of the coldest x voxels.
-					@return Return ose value in Gy.
-				*/
-				DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute) const;
-			};
-
-		}
+		{
+		public:
+			enum complexStatistics { Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx};
+			typedef boost::shared_ptr<std::vector<std::pair<DoseTypeGy, VoxelGridID> > > ResultListPointer;
+			typedef boost::shared_ptr<DoseStatistics> DoseStatisticsPointer;
+			typedef std::map<DoseTypeGy, VolumeType> DoseToVolumeFunctionType;
+			typedef std::map<VolumeType, DoseTypeGy> VolumeToDoseFunctionType;
+
+		private:
+			double getValue(std::map<double, double> aMap, double key, bool findNearestValueInstead,
+			                double& storedKey) const;
+
+			std::map<double, double>::const_iterator findNearestKeyInMap(std::map<double, double> aMap, double key) const;
+
+			DoseStatisticType _maximum;
+			DoseStatisticType _minimum;
+			ResultListPointer _maximumVoxelPositions;
+			ResultListPointer _minimumVoxelPositions;
+			DoseStatisticType _mean;
+			DoseStatisticType _stdDeviation;
+			unsigned int _numVoxels;
+			VolumeType _volume;
+			DoseToVolumeFunctionType _Dx;
+			VolumeToDoseFunctionType _Vx;
+			VolumeToDoseFunctionType _MOHx;
+			VolumeToDoseFunctionType _MOCx;
+			VolumeToDoseFunctionType _MaxOHx;
+			VolumeToDoseFunctionType _MinOCx;
+
+		public:
+			/*! @brief Standard Constructor
+			*/
+			//DoseStatistics();
+
+			/*! @brief Constructor
+				@detail the dose statistic values are set. Complex values maximumVoxelLocation, maximumVoxelLocation, Dx, Vx, MOHx, MOCx, MaxOHx and MinOCx are optional
+			*/
+			DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean,
+			               DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume,
+			               ResultListPointer minimumVoxelPositions = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >
+			               (std::vector<std::pair<DoseTypeGy, VoxelGridID> >()),
+			               ResultListPointer maximumVoxelPositions = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >
+			               (std::vector<std::pair<DoseTypeGy, VoxelGridID> >()),
+			               DoseToVolumeFunctionType Dx = DoseToVolumeFunctionType(),
+			               VolumeToDoseFunctionType Vx = VolumeToDoseFunctionType(),
+			               VolumeToDoseFunctionType MOHx = VolumeToDoseFunctionType(),
+			               VolumeToDoseFunctionType MOCx = VolumeToDoseFunctionType(),
+			               VolumeToDoseFunctionType MaxOHx = VolumeToDoseFunctionType(),
+			               VolumeToDoseFunctionType MinOCx = VolumeToDoseFunctionType());
+
+			~DoseStatistics();
+
+			void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions);
+			void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions);
+			void setDx(const DoseToVolumeFunctionType& DxValues);
+			void setVx(const VolumeToDoseFunctionType& VxValues);
+			void setMOHx(const VolumeToDoseFunctionType& MOHxValues);
+			void setMOCx(const VolumeToDoseFunctionType& MOCxValues);
+			void setMaxOHx(const VolumeToDoseFunctionType& MaxOHxValues);
+			void setMinOCx(const VolumeToDoseFunctionType& MinOCxValues);
+
+			/*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy.
+			*/
+			unsigned int getNumberOfVoxels() const;
+
+			VolumeType getVolume() const;
+
+			/*! @brief Get the maximum of the current dose distribution.
+				@return Return the maximum dose in Gy
+			*/
+			DoseStatisticType getMaximum() const;
+
+			/*! @brief Get a vector of the the maximum dose VoxelGridIDs together with their dose value in Gy
+				@exception InvalidDoseException if the vector has not been set (i.e. is empty)
+			*/
+			ResultListPointer getMaximumPositions() const;
+
+			/*! @brief Get the minimum of the current dose distribution.
+				@return Return the minimum dose in Gy
+			*/
+			DoseStatisticType getMinimum() const;
+
+			/*! @brief Get a vector of the the minimum dose VoxelGridIDs together with their dose value in Gy
+			@exception InvalidDoseException if the vector has not been set (i.e. is empty)
+			*/
+			ResultListPointer getMinimumPositions() const;
+
+			/*! @brief Get the mean of the current dose distribution.
+				@return Return the mean dose in Gy
+			*/
+			DoseStatisticType getMean() const;
+
+			/*! @brief Get the standard deviation of the current dose distribution.
+				@return Return the standard deviation in Gy
+			*/
+			DoseStatisticType getStdDeviation() const;
+
+			/*! @brief Get the variance of of the current dose distribution.
+				@return Return the variance in Gy
+			*/
+			DoseStatisticType getVariance() const;
+
+			/*! @brief Get Vx: the volume irradiated with a dose >= x.
+				@return Return absolute volume in absolute cm^3.
+				@exception NoDataException if the Vx values have not been set (i.e. the vector is empty)
+				@exception NoDataException if the requested Dose is not in the vector
+			*/
+			VolumeType getVx(DoseTypeGy xDoseAbsolute) const;
+			VolumeType getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue,
+			                 DoseTypeGy& nearestXDose) const;
+			DoseToVolumeFunctionType getAllVx() const;
+
+			/*! @brief Get Dx: the minimal dose delivered to part x of the current volume.
+				@return Return dose value in Gy.
+				@exception InvalidDoseException if the Dx values have not been set (i.e. the vector is empty)
+			*/
+			DoseTypeGy getDx(VolumeType xVolumeAbsolute, bool findNearestValue,
+			                 VolumeType& nearestXVolume) const;
+			DoseTypeGy getDx(VolumeType xVolumeAbsolute) const;
+			VolumeToDoseFunctionType getAllDx() const;
+
+			/*! @brief Get MOHx: mean dose of the hottest x voxels.
+				@return Return dose value in Gy.
+				@exception InvalidDoseException if the values have not been set (i.e. the vector is empty)
+			*/
+			DoseTypeGy getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue,
+			                   VolumeType& nearestXVolume) const;
+			DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) const;
+			VolumeToDoseFunctionType getAllMOHx() const;
+
+			/*! @brief Get MOCx: mean dose of the coldest x voxels.
+				@return Return dose value in Gy.
+				@exception InvalidDoseException if the values have not been set (i.e. the vector is empty)
+			*/
+			DoseTypeGy getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue,
+			                   VolumeType& nearestXVolume) const;
+			DoseTypeGy getMOCx(VolumeType xVolumeAbsolute) const;
+			VolumeToDoseFunctionType getAllMOCx() const;
+
+			/*! @brief Get MaxOHx: Maximum outside of the hottest x voxels.
+				@return Return dose value in Gy.
+				@exception InvalidDoseException if the values have not been set (i.e. the vector is empty)
+			*/
+			DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue,
+			                     VolumeType& nearestXVolume) const;
+			DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute) const;
+			VolumeToDoseFunctionType getAllMaxOHx() const;
+
+			/*! @brief Get MinOCx: Minimum outside of the coldest x voxels.
+				@return Return dose value in Gy.
+				@exception InvalidDoseException if the values have not been set (i.e. the vector is empty)
+			*/
+			DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue,
+			                     VolumeType& nearestXVolume) const;
+			DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute) const;
+			VolumeToDoseFunctionType getAllMinOCx() const;
+		};
+
 	}
+}
 
 
 #endif
diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp
new file mode 100644
index 0000000..2039fcf
--- /dev/null
+++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp
@@ -0,0 +1,554 @@
+// -----------------------------------------------------------------------
+// RTToolbox - DKFZ radiotherapy quantitative evaluation library
+//
+// Copyright (c) German Cancer Research Center (DKFZ),
+// Software development for Integrated Diagnostics and Therapy (SIDT).
+// ALL RIGHTS RESERVED.
+// See rttbCopyright.txt or
+// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
+//
+// This software is distributed WITHOUT ANY WARRANTY; without even
+// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+// PURPOSE.  See the above copyright notices for more information.
+//
+//------------------------------------------------------------------------
+/*!
+// @file
+// @version $Revision: 877 $ (last changed revision)
+// @date    $Date: 2015-01-09 10:51:10 +0100 (Fr, 09 Jan 2015) $ (last change date)
+// @author  $Author: hentsch $ (last changed by)
+*/
+
+#include "rttbDoseStatisticsCalculator.h"
+
+#include <boost/shared_ptr.hpp>
+#include <boost/make_shared.hpp>
+#include <boost/assign/list_of.hpp>
+
+#include "rttbNullPointerException.h"
+#include "rttbInvalidDoseException.h"
+#include "rttbInvalidParameterException.h"
+
+namespace rttb
+{
+
+	namespace algorithms
+	{
+		DoseStatisticsCalculator::DoseStatisticsCalculator()
+		{
+			simpleDoseStatisticsCalculated = false;
+			_doseIterator = NULL;
+		}
+
+		DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator)
+		{
+			setDoseIterator(aDoseIterator);
+			simpleDoseStatisticsCalculated = false;
+		}
+
+
+		DoseStatisticsCalculator::~DoseStatisticsCalculator()
+		{
+
+		}
+
+
+		void DoseStatisticsCalculator::setDoseIterator(DoseIteratorPointer aDoseIterator)
+		{
+			if (aDoseIterator == NULL)
+			{
+				throw core::NullPointerException("DoseIterator must not be NULL");
+			}
+			else
+			{
+				_doseIterator = aDoseIterator;
+			}
+		}
+
+		DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const
+		{
+			return _doseIterator;
+		}
+
+
+		DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics(
+		    bool computeComplexMeasures, const std::vector<double>& precomputeDoseValues,
+		    const std::vector<double>& precomputeVolumeValues)
+		{
+
+			if (!_doseIterator)
+			{
+				throw core::NullPointerException("_doseIterator must not be NULL!");
+			}
+
+			//"simple" dose statistics are obligatory
+			calculateSimpleDoseStatistics();
+
+			if (computeComplexMeasures)
+			{
+				//more complex dose statistics are optional
+				calculateComplexDoseStatistics(precomputeDoseValues, precomputeVolumeValues);
+			}
+
+			return _statistics;
+
+		}
+
+
+		void DoseStatisticsCalculator::calculateSimpleDoseStatistics()
+		{
+			_doseVector.clear();
+			_voxelProportionVector.clear();
+
+			std::multimap<double, int> doseValueVSIndexMap;
+			std::vector<double> voxelProportionVectorTemp;
+
+
+			DoseStatisticType maximumDose = 0;
+			DoseStatisticType minimumDose = 0;
+			DoseStatisticType meanDose = 0;
+			DoseStatisticType stdDeviationDose = 0;
+
+
+			float sum = 0;
+			unsigned int numVoxels = 0;
+			float squareSum = 0;
+			VolumeType volume;
+
+			_doseIterator->reset();
+			int i = 0;
+			DoseTypeGy doseValue = 0;
+
+			while (_doseIterator->isPositionValid())
+			{
+				doseValue = _doseIterator->getCurrentDoseValue();
+
+				if (i == 0)
+				{
+					minimumDose = doseValue;
+					volume = _doseIterator->getCurrentVoxelVolume();
+				}
+
+				rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction();
+				sum += doseValue * voxelProportion;
+				numVoxels += voxelProportion;
+				squareSum += doseValue * doseValue * voxelProportion;
+
+				if (doseValue > maximumDose)
+				{
+					maximumDose = doseValue;
+				}
+				else if (doseValue < minimumDose)
+				{
+					minimumDose = doseValue;
+				}
+
+				voxelProportionVectorTemp.push_back(voxelProportion);
+				doseValueVSIndexMap.insert(std::pair<double, int>(doseValue, i));
+
+				i++;
+				_doseIterator->next();
+			}
+
+			if (numVoxels != 0)
+			{
+				meanDose = sum / numVoxels;
+				DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose);
+
+				//standard deviation is defined only for n>=2
+				if (varianceDose < errorConstant || numVoxels <= 2)
+				{
+					stdDeviationDose = 0;
+				}
+				else
+				{
+					stdDeviationDose = pow(varianceDose, 0.5);
+				}
+
+			}
+
+			//sort dose values and corresponding volume fractions in member variables
+			std::multimap<double, int>::iterator it;
+
+			for (it = doseValueVSIndexMap.begin(); it != doseValueVSIndexMap.end(); ++it)
+			{
+				_doseVector.push_back((float)(*it).first);
+				_voxelProportionVector.push_back(voxelProportionVectorTemp.at((*it).second));
+			}
+
+			volume *= numVoxels;
+
+			_statistics = boost::make_shared<DoseStatistics>(minimumDose, maximumDose, meanDose, stdDeviationDose, numVoxels,
+			              volume);
+
+			simpleDoseStatisticsCalculated = true;
+
+			ResultListPointer minimumVoxelPositions = computeMinimumPositions(100);
+			ResultListPointer maximumVoxelPositions = computeMaximumPositions(100);
+
+			_statistics->setMinimumVoxelPositions(minimumVoxelPositions);
+			_statistics->setMaximumVoxelPositions(maximumVoxelPositions);
+
+
+		}
+
+
+		void DoseStatisticsCalculator::calculateComplexDoseStatistics(const std::vector<double>& precomputeDoseValues,
+		        const std::vector<double>& precomputeVolumeValues)
+		{
+			std::vector<double> precomputeDoseValuesNonConst = precomputeDoseValues;
+			std::vector<double> precomputeVolumeValuesNonConst = precomputeVolumeValues;
+
+			//set default values
+			if (precomputeDoseValues.empty())
+			{
+				std::vector<double> defaultPrecomputeDoseValues = boost::assign::list_of(0.02 * _statistics->getMaximum())(
+				            0.05 * _statistics->getMaximum())(0.1 * _statistics->getMaximum())(0.9 * _statistics->getMaximum())(
+				            0.95 * _statistics->getMaximum())(0.98 * _statistics->getMaximum());
+				precomputeDoseValuesNonConst = defaultPrecomputeDoseValues;
+			}
+
+			if (precomputeVolumeValues.empty())
+			{
+				std::vector<double> defaultPrecomputeVolumeValues = boost::assign::list_of(0.02 * _statistics->getNumberOfVoxels())(
+				            0.05 * _statistics->getNumberOfVoxels())(0.1 * _statistics->getNumberOfVoxels())(0.9 * _statistics->getNumberOfVoxels())
+				        (0.95 * _statistics->getNumberOfVoxels())(0.98 * _statistics->getNumberOfVoxels());
+				precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues;
+			}
+
+			DoseToVolumeFunctionType Vx = computeDoseToVolumeMulti(precomputeDoseValuesNonConst, DoseStatistics::Vx);
+			VolumeToDoseFunctionType Dx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::Dx);
+			VolumeToDoseFunctionType MOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOHx);
+			VolumeToDoseFunctionType MOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOCx);
+			VolumeToDoseFunctionType MaxOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst,
+			                                  DoseStatistics::MaxOHx);
+			VolumeToDoseFunctionType MinOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst,
+			                                  DoseStatistics::MinOCx);
+
+			_statistics->setVx(Vx);
+			_statistics->setDx(Dx);
+			_statistics->setMOHx(MOHx);
+			_statistics->setMOCx(MOCx);
+			_statistics->setMaxOHx(MaxOHx);
+			_statistics->setMinOCx(MinOCx);
+		}
+
+
+
+		DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions(
+		    unsigned int maxNumberMaxima)
+		{
+			if (!simpleDoseStatisticsCalculated)
+			{
+				throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumAndPosition()");
+
+			}
+
+			ResultListPointer maxVoxelVector = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+
+			unsigned int count = 0;
+			this->_doseIterator->reset();
+			DoseTypeGy doseValue = 0;
+
+			while (_doseIterator->isPositionValid() && count < maxNumberMaxima)
+			{
+				doseValue = _doseIterator->getCurrentDoseValue();
+
+				if (doseValue == _statistics->getMaximum())
+				{
+					VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID();
+					std::pair<DoseTypeGy, VoxelGridID> voxel(doseValue, currentID);
+					maxVoxelVector->push_back(voxel);
+					count++;
+				}
+
+				_doseIterator->next();
+			}
+
+			return maxVoxelVector;
+		}
+
+		DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMinimumPositions(
+		    unsigned int maxNumberMinima)
+		{
+			if (!simpleDoseStatisticsCalculated)
+			{
+				throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumAndPosition()");
+
+			}
+
+			ResultListPointer minVoxelVector = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+
+			/*! @todo: Architecture Annotation:
+				Finding the positions for the minimum only once reduces computation time,
+				but will require sensible use by the programmers. To be save the output vector
+				minVoxelVector will be always cleared here to garantee that no false values are
+				presented. This change may be revoced to increase computation speed later on
+				(only compute if(minVoxelVector->size()==0)).
+			*/
+			unsigned int count = 0;
+			this->_doseIterator->reset();
+			DoseTypeGy doseValue = 0;
+
+			while (_doseIterator->isPositionValid() && count < maxNumberMinima)
+			{
+				doseValue = _doseIterator->getCurrentDoseValue();
+
+				if (doseValue == _statistics->getMinimum())
+				{
+					VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID();
+					std::pair<DoseTypeGy, VoxelGridID> voxel(doseValue, currentID);
+					minVoxelVector->push_back(voxel);
+					count++;
+				}
+
+				_doseIterator->next();
+			}
+
+			return minVoxelVector;
+		}
+
+		VolumeType DoseStatisticsCalculator::computeVx(DoseTypeGy xDoseAbsolute) const
+		{
+			rttb::FractionType count = 0;
+			_doseIterator->reset();
+
+			DoseTypeGy currentDose = 0;
+
+			while (_doseIterator->isPositionValid())
+			{
+				currentDose = _doseIterator->getCurrentDoseValue();
+
+				if (currentDose >= xDoseAbsolute)
+				{
+					count += _doseIterator->getCurrentRelevantVolumeFraction();
+				}
+
+				_doseIterator->next();
+			}
+
+			return count * this->_doseIterator->getCurrentVoxelVolume();
+		}
+
+		DoseTypeGy DoseStatisticsCalculator::computeDx(VolumeType xVolumeAbsolute) const
+		{
+			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+			DoseTypeGy resultDose = 0;
+
+			double countVoxels = 0;
+			int i = _doseVector.size() - 1;
+
+			for (; i >= 0; i--)
+			{
+				countVoxels += _voxelProportionVector.at(i);
+
+				if (countVoxels >= noOfVoxel)
+				{
+					break;
+				}
+			}
+
+			if (i >= 0)
+			{
+				resultDose = _doseVector.at(i);
+			}
+			else
+			{
+				resultDose = _statistics->getMinimum();
+			}
+
+			return resultDose;
+		}
+
+		DoseTypeGy DoseStatisticsCalculator::computeMOHx(VolumeType xVolumeAbsolute) const
+		{
+			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+
+
+			if (noOfVoxel == 0)
+			{
+				return 0;
+			}
+			else
+			{
+				double countVoxels = 0;
+				double sum = 0;
+
+				for (int i = _doseVector.size() - 1; i >= 0; i--)
+				{
+					double voxelProportion = _voxelProportionVector.at(i);
+					countVoxels += voxelProportion;
+					sum += _doseVector.at(i) * voxelProportion;
+
+					if (countVoxels >= noOfVoxel)
+					{
+						break;
+					}
+				}
+
+				return (DoseTypeGy)(sum / noOfVoxel);
+			}
+		}
+
+		DoseTypeGy DoseStatisticsCalculator::computeMOCx(DoseTypeGy xVolumeAbsolute) const
+		{
+			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+
+
+			if (noOfVoxel == 0)
+			{
+				return 0;
+			}
+			else
+			{
+				double countVoxels = 0;
+				double sum = 0;
+				std::vector<DoseTypeGy>::const_iterator it = _doseVector.begin();
+				std::vector<double>::const_iterator itD = _voxelProportionVector.begin();
+
+				for (; it != _doseVector.end(); ++it, ++itD)
+				{
+					double voxelProportion = *itD;
+					countVoxels += voxelProportion;
+					sum += (*it) * voxelProportion;
+
+					if (countVoxels >= noOfVoxel)
+					{
+						break;
+					}
+				}
+
+				return (DoseTypeGy)(sum / noOfVoxel);
+			}
+		}
+
+		DoseTypeGy DoseStatisticsCalculator::computeMaxOHx(DoseTypeGy xVolumeAbsolute) const
+		{
+			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+			DoseTypeGy resultDose = 0;
+
+			double countVoxels = 0;
+			int i = _doseVector.size() - 1;
+
+			for (; i >= 0; i--)
+			{
+				countVoxels += _voxelProportionVector.at(i);
+
+				if (countVoxels >= noOfVoxel)
+				{
+					break;
+				}
+			}
+
+			if (i - 1 >= 0)
+			{
+				resultDose = _doseVector.at(i - 1);
+			}
+
+			return resultDose;
+		}
+
+		DoseTypeGy DoseStatisticsCalculator::computeMinOCx(DoseTypeGy xVolumeAbsolute) const
+		{
+			double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume();
+			DoseTypeGy resultDose = 0;
+
+			double countVoxels = 0;
+			std::vector<DoseTypeGy>::const_iterator it = _doseVector.begin();
+			std::vector<double>::const_iterator itD = _voxelProportionVector.begin();
+
+			for (; itD != _voxelProportionVector.end(); ++itD, ++it)
+			{
+				countVoxels += *itD;
+
+				if (countVoxels >= noOfVoxel)
+				{
+					break;
+				}
+			}
+
+			if (it != _doseVector.end())
+			{
+				++it;
+
+				if (it != _doseVector.end())
+				{
+					resultDose = *it;
+				}
+				else
+				{
+					resultDose = (DoseTypeGy)_statistics->getMaximum();
+				}
+			}
+			else
+			{
+				resultDose = (DoseTypeGy)_statistics->getMinimum();
+			}
+
+			return resultDose;
+		}
+
+		DoseStatisticsCalculator::DoseToVolumeFunctionType DoseStatisticsCalculator::computeDoseToVolumeMulti(
+		    const std::vector<double>& precomputeDoseValues, DoseStatistics::complexStatistics name) const
+		{
+			DoseToVolumeFunctionType VxMulti;
+			DoseTypeGy maxDose = _statistics->getMaximum();
+
+			for (int i = 0; i < precomputeDoseValues.size(); ++i)
+			{
+				if (name == DoseStatistics::Vx)
+				{
+					VxMulti.insert(std::pair<DoseTypeGy, VolumeType>(precomputeDoseValues.at(i),
+					               computeVx(maxDose * (precomputeDoseValues.at(i) / 100.0))));
+				}
+			}
+
+			return VxMulti;
+		}
+
+		DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti(
+		    const std::vector<double>& precomputeVolumeValues, DoseStatistics::complexStatistics name) const
+		{
+			VolumeToDoseFunctionType multiValues;
+			VolumeType maxVolume = _statistics->getNumberOfVoxels();
+
+			for (int i = 0; i < precomputeVolumeValues.size(); ++i)
+			{
+				switch (name)
+				{
+					case DoseStatistics::Dx:
+						multiValues.insert(std::pair<VolumeType, DoseTypeGy>(precomputeVolumeValues.at(i),
+						                   computeDx(maxVolume * (precomputeVolumeValues.at(i) / 100.0))));
+						break;
+
+					case DoseStatistics::MOHx:
+						multiValues.insert(std::pair<VolumeType, DoseTypeGy>(precomputeVolumeValues.at(i),
+						                   computeMOHx(maxVolume * (precomputeVolumeValues.at(i) / 100.0))));
+						break;
+
+					case DoseStatistics::MOCx:
+						multiValues.insert(std::pair<VolumeType, DoseTypeGy>(precomputeVolumeValues.at(i),
+						                   computeMOCx(maxVolume * (precomputeVolumeValues.at(i) / 100.0))));
+						break;
+
+					case DoseStatistics::MaxOHx:
+						multiValues.insert(std::pair<VolumeType, DoseTypeGy>(precomputeVolumeValues.at(i),
+						                   computeMaxOHx(maxVolume * (precomputeVolumeValues.at(i) / 100.0))));
+						break;
+
+					case DoseStatistics::MinOCx:
+						multiValues.insert(std::pair<VolumeType, DoseTypeGy>(precomputeVolumeValues.at(i),
+						                   computeMinOCx(maxVolume * (precomputeVolumeValues.at(i) / 100.0))));
+						break;
+
+					default:
+						throw core::InvalidParameterException("unknown DoseStatistics name!");
+				}
+			}
+
+			return multiValues;
+		}
+
+	}//end namespace algorithms
+}//end namespace rttb
+
diff --git a/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h
new file mode 100644
index 0000000..8d01d25
--- /dev/null
+++ b/code/algorithms/rttbDoseStatisticsCalculator.h
@@ -0,0 +1,120 @@
+// -----------------------------------------------------------------------
+// RTToolbox - DKFZ radiotherapy quantitative evaluation library
+//
+// Copyright (c) German Cancer Research Center (DKFZ),
+// Software development for Integrated Diagnostics and Therapy (SIDT).
+// ALL RIGHTS RESERVED.
+// See rttbCopyright.txt or
+// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
+//
+// This software is distributed WITHOUT ANY WARRANTY; without even
+// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+// PURPOSE.  See the above copyright notices for more information.
+//
+//------------------------------------------------------------------------
+/*!
+// @file
+// @version $Revision: 707 $ (last changed revision)
+// @date    $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date)
+// @author  $Author: floca $ (last changed by)
+*/
+
+#ifndef __DOSE_STATISTICS_CALCULATOR_H
+#define __DOSE_STATISTICS_CALCULATOR_H
+
+#include <vector>
+
+#include "rttbBaseType.h"
+#include "rttbDoseIteratorInterface.h"
+#include "rttbDoseStatistics.h"
+
+namespace rttb
+{
+
+	namespace algorithms
+	{
+
+		/*! @class DoseStatisticsCalculator
+		@brief This is a class calculating different statistical values from a rt dose distribution
+		These values range from standard statistical values such as minimum, maximum and mean to more
+		dose specific properties such as Vx (volume irradiated with a dose >=x), Dx (minimal dose delivered
+		to x% of the VOI, and MOHx (mean in the hottest volume).
+		*/
+		class DoseStatisticsCalculator
+		{
+		public:
+			typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer;
+			typedef DoseStatistics::ResultListPointer ResultListPointer;
+			typedef DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer;
+			typedef DoseStatistics::DoseToVolumeFunctionType DoseToVolumeFunctionType;
+			typedef DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType;
+
+		private:
+			DoseIteratorPointer _doseIterator;
+
+			/* Contains relevant dose values sorted in descending order.
+			*/
+			std::vector<DoseTypeGy> _doseVector;
+
+			/*! Contains the corresponding voxel proportions to the values in doseVector.
+			*/
+			std::vector<double> _voxelProportionVector;
+
+
+			DoseStatisticsPointer _statistics;
+
+			bool simpleDoseStatisticsCalculated;
+
+			/*! @brief Calculation of basic dose statistics. It will be called in Constructor.
+				@exception NullPointerException Thrown if _doseIterator is NULL.
+			*/
+
+
+			ResultListPointer computeMaximumPositions(unsigned int maxNumberMaxima);
+
+			ResultListPointer computeMinimumPositions(unsigned int maxNumberMinima);
+			VolumeType computeVx(DoseTypeGy xDoseAbsolute) const;
+			DoseTypeGy computeDx(VolumeType xVolumeAbsolute) const;
+			DoseTypeGy computeMOHx(VolumeType xVolumeAbsolute) const;
+			DoseTypeGy computeMOCx(DoseTypeGy xVolumeAbsolute) const;
+			DoseTypeGy computeMaxOHx(DoseTypeGy xVolumeAbsolute) const;
+			DoseTypeGy computeMinOCx(DoseTypeGy xVolumeAbsolute) const;
+
+			DoseToVolumeFunctionType computeDoseToVolumeMulti(const std::vector<double>& precomputeDoseValues,
+			        DoseStatistics::complexStatistics name) const;
+			VolumeToDoseFunctionType computeVolumeToDoseFunctionMulti(const std::vector<double>& precomputeVolumeValues,
+			        DoseStatistics::complexStatistics name) const;
+
+			void calculateSimpleDoseStatistics();
+			void calculateComplexDoseStatistics(const std::vector<double>& precomputeDoseValues,
+			                                    const std::vector<double>& precomputeVolumeValues);
+
+
+		public:
+			/*! @brief Standard Constructor
+			*/
+			DoseStatisticsCalculator();
+
+			~DoseStatisticsCalculator();
+
+			/*! @brief Constructor
+			*/
+			DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator);
+
+			/*! @brief Set new dose data for statistics. Statistics will be re-initialized.
+			*/
+			void setDoseIterator(DoseIteratorPointer aDoseIterator);
+
+			DoseIteratorPointer getDoseIterator() const;
+
+			DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false,
+			        const std::vector<double>& precomputeDoseValues = std::vector<double>(),
+			        const std::vector<double>& precomputeVolumeValues = std::vector<double>());
+
+		};
+
+	}
+}
+
+
+#endif
diff --git a/code/core/files.cmake b/code/core/files.cmake
index 069db8c..56c703a 100644
--- a/code/core/files.cmake
+++ b/code/core/files.cmake
@@ -1,58 +1,60 @@
 SET(CPP_FILES 
   rttbDoseIteratorInterface.cpp
   rttbDVH.cpp
   rttbDVHCalculator.cpp
   rttbDVHSet.cpp
+  rttbDataNotAvailableException.cpp
   rttbException.cpp		
   rttbGenericDoseIterator.cpp
   rttbGenericMaskedDoseIterator.cpp
   rttbGeometricInfo.cpp
   rttbIndexOutOfBoundsException.cpp
   rttbInvalidDoseException.cpp
   rttbInvalidParameterException.cpp
   rttbMappingOutsideOfImageException.cpp
   rttbMaskedDoseIteratorInterface.cpp
   rttbMaskVoxel.cpp
   rttbNullPointerException.cpp
   rttbPaddingException.cpp
   rttbPhysicalInfo.cpp
   rttbStructure.cpp
   rttbStructureSet.cpp
   rttbStrVectorStructureSetGenerator.cpp
   )
 
 SET(H_FILES 
   rttbBaseType.h
+  rttbDataNotAvailableException.h
   rttbDoseAccessorGeneratorBase.h
   rttbDoseAccessorGeneratorInterface.h
   rttbDoseAccessorInterface.h
   rttbDoseIteratorInterface.h
   rttbDVH.h
   rttbDVHCalculator.h
   rttbDVHGeneratorInterface.h
   rttbDVHSet.h
   rttbException.h
   rttbExceptionMacros.h
   rttbGenericDoseIterator.h
   rttbGenericMaskedDoseIterator.h
   rttbGeometricInfo.h
   rttbIndexConversionInterface.h
   rttbIndexOutOfBoundsException.h
   rttbInvalidDoseException.h
   rttbInvalidParameterException.h
   rttbMappingOutsideOfImageException.h
   rttbMaskAccessorGeneratorBase.h
   rttbMaskAccessorGeneratorInterface.h
   rttbMaskAccessorInterface.h
   rttbMaskedDoseIteratorInterface.h
   rttbMaskVoxel.h
   rttbMutableDoseAccessorInterface.h
   rttbMutableMaskAccessorInterface.h
   rttbNullPointerException.h
   rttbPaddingException.h
   rttbPhysicalInfo.h
   rttbStructure.h
   rttbStructureSet.h
   rttbStructureSetGeneratorInterface.h
   rttbStrVectorStructureSetGenerator.h
 )
diff --git a/code/core/rttbDataNotAvailableException.cpp b/code/core/rttbDataNotAvailableException.cpp
new file mode 100644
index 0000000..357c781
--- /dev/null
+++ b/code/core/rttbDataNotAvailableException.cpp
@@ -0,0 +1,43 @@
+// -----------------------------------------------------------------------
+// RTToolbox - DKFZ radiotherapy quantitative evaluation library
+//
+// Copyright (c) German Cancer Research Center (DKFZ),
+// Software development for Integrated Diagnostics and Therapy (SIDT).
+// ALL RIGHTS RESERVED.
+// See rttbCopyright.txt or
+// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
+//
+// This software is distributed WITHOUT ANY WARRANTY; without even
+// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+// PURPOSE.  See the above copyright notices for more information.
+//
+//------------------------------------------------------------------------
+/*!
+// @file
+// @version $Revision: 741 $ (last changed revision)
+// @date    $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date)
+// @author  $Author: hentsch $ (last changed by)
+*/
+
+#include <iostream>
+#include <sstream>
+
+#include "rttbDataNotAvailableException.h"
+
+namespace rttb
+{
+	namespace core
+	{
+
+		const char* DataNotAvailableException::what() const throw()
+		{
+			return rttb_what.c_str();
+		}
+
+		const char* DataNotAvailableException::GetNameOfClass() const
+		{
+			return "DataNotAvailableException";
+		}
+
+	}//end namespace core
+}//end namespace rttb
diff --git a/code/core/rttbDataNotAvailableException.h b/code/core/rttbDataNotAvailableException.h
new file mode 100644
index 0000000..bfa78c1
--- /dev/null
+++ b/code/core/rttbDataNotAvailableException.h
@@ -0,0 +1,58 @@
+// -----------------------------------------------------------------------
+// RTToolbox - DKFZ radiotherapy quantitative evaluation library
+//
+// Copyright (c) German Cancer Research Center (DKFZ),
+// Software development for Integrated Diagnostics and Therapy (SIDT).
+// ALL RIGHTS RESERVED.
+// See rttbCopyright.txt or
+// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
+//
+// This software is distributed WITHOUT ANY WARRANTY; without even
+// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+// PURPOSE.  See the above copyright notices for more information.
+//
+//------------------------------------------------------------------------
+/*!
+// @file
+// @version $Revision: 741 $ (last changed revision)
+// @date    $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date)
+// @author  $Author: hentsch $ (last changed by)
+*/
+
+#ifndef __DATA_NOT_AVAILABLE_EXCEPTION_H
+#define __DATA_NOT_AVAILABLE_EXCEPTION_H
+
+#include <string>
+#include <exception>
+
+#include "rttbException.h"
+
+
+namespace rttb
+{
+	namespace core
+	{
+
+		/*! @class DataNotAvailableException
+			@brief This exception will be thrown if the requested data is not available.
+		*/
+		class DataNotAvailableException : public Exception
+		{
+		public:
+			DataNotAvailableException(const std::string& aWhat) : Exception(aWhat) {}
+
+			virtual ~DataNotAvailableException() throw() {}
+
+			/*! @brief Get the exception description
+			*/
+			virtual const char* what() const throw();
+
+			/*! @brief Get the name of the exception class
+			*/
+			virtual const char* GetNameOfClass() const;
+		};
+
+	}
+}
+
+#endif
diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp
index 88beb4e..39bf48e 100644
--- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp
+++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp
@@ -1,380 +1,464 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE.  See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date    $Date$ (last change date)
 // @author  $Author$ (last changed by)
 */
 
 #include <ostream>
 #include <fstream>
 
 #include "rttbDoseStatisticsXMLWriter.h"
 #include "rttbNullPointerException.h"
 #include "rttbInvalidParameterException.h"
+#include "rttbDataNotAvailableException.h"
 #include "rttbException.h"
 
 #include "boost/lexical_cast.hpp"
 
 namespace rttb
 {
 	namespace io
 	{
 		namespace other
 		{
 			static const std::string xmlattrTag = "<xmlattr>.x";
 			static const std::string statisticsTag = "statistics";
 			static const std::string columnSeparator = "@";
 
 			boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics,
 			        DoseTypeGy aReferenceDose)
 			{
 
 				using boost::property_tree::ptree;
 				ptree pt;
 
-				rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator =
-				    aDoseStatistics->getDoseIterator();
-				rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor =
-				    spDoseIterator->getDoseAccessor();
+				/*rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator =
+					aDoseStatistics->getDoseIterator();
+					rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor =
+					spDoseIterator->getDoseAccessor();
 
-				boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs =
-				    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
-				ResultListPointer spMyResultPairs(myResultPairs);
-				boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs2 =
-				    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
-				ResultListPointer spMyResultPairs2(myResultPairs2);
+					boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs =
+					boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+					ResultListPointer spMyResultPairs(myResultPairs);
+					boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs2 =
+					boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+					ResultListPointer spMyResultPairs2(myResultPairs2);*/
 
 				pt.put(statisticsTag + ".numberOfVoxels", aDoseStatistics->getNumberOfVoxels());
-				pt.put(statisticsTag + ".volume", aDoseStatistics->getVx(0));
-				pt.put(statisticsTag + ".minimum", aDoseStatistics->getMinimum(spMyResultPairs));
-				std::vector<std::pair<DoseTypeGy, VoxelGridID> >::iterator pairIt = spMyResultPairs->begin();
+				pt.put(statisticsTag + ".volume", aDoseStatistics->getVolume());
+				pt.put(statisticsTag + ".minimum", aDoseStatistics->getMinimum());
+
+				auto minimumPositions = aDoseStatistics->getMinimumPositions();
+				std::vector<std::pair<DoseTypeGy, VoxelGridID> >::iterator pairItMin = minimumPositions->begin();
 				int count = 0;
 
-				for (; pairIt != spMyResultPairs->end() && count < 100 ; ++pairIt) //output max. 100 minimum
+				for (; pairItMin != minimumPositions->end() && count < 100; ++pairItMin) //output max. 100 minimum
 				{
-					VoxelGridID voxelID = pairIt->second;
+					VoxelGridID voxelID = pairItMin->second;
 					pt.add(statisticsTag + ".minimum.voxelGridID", voxelID);
 
-					VoxelGridIndex3D voxelIndex3D;
-					spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D);
-					std::string voxelIndex3DStr;
-					voxelIndex3DStr = boost::lexical_cast<std::string>(voxelIndex3D.x()) + ","
-					                  + boost::lexical_cast<std::string>(voxelIndex3D.y()) + ","
-					                  + boost::lexical_cast<std::string>(voxelIndex3D.z());
-					pt.add(statisticsTag + ".minimum.voxelIndex3D", voxelIndex3DStr);
-
-					WorldCoordinate3D worldCoor;
-					spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor);
-					std::string worldCoorStr;
-					worldCoorStr = boost::lexical_cast<std::string>(worldCoor.x()) + ","
-					               + boost::lexical_cast<std::string>(worldCoor.y()) + ","
-					               + boost::lexical_cast<std::string>(worldCoor.z());
-					pt.add(statisticsTag + ".minimum.worldCoordinate", worldCoorStr);
+					//VoxelGridIndex3D voxelIndex3D;
+					//spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D);
+					//std::string voxelIndex3DStr;
+					//voxelIndex3DStr = boost::lexical_cast<std::string>(voxelIndex3D.x()) + ","
+					//                  + boost::lexical_cast<std::string>(voxelIndex3D.y()) + ","
+					//                  + boost::lexical_cast<std::string>(voxelIndex3D.z());
+					//pt.add(statisticsTag + ".minimum.voxelIndex3D", voxelIndex3DStr);
+
+					//WorldCoordinate3D worldCoor;
+					//spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor);
+					//std::string worldCoorStr;
+					//worldCoorStr = boost::lexical_cast<std::string>(worldCoor.x()) + ","
+					//               + boost::lexical_cast<std::string>(worldCoor.y()) + ","
+					//               + boost::lexical_cast<std::string>(worldCoor.z());
+					//pt.add(statisticsTag + ".minimum.worldCoordinate", worldCoorStr);
 
 					count++;
 				}
 
-				pt.put(statisticsTag + ".maximum", aDoseStatistics->getMaximum(spMyResultPairs2));
-				std::vector<std::pair<DoseTypeGy, VoxelGridID> >::iterator pairIt2 = spMyResultPairs2->begin();
+				pt.put(statisticsTag + ".maximum", aDoseStatistics->getMaximum());
+
+				auto maximumPositions = aDoseStatistics->getMaximumPositions();
+				std::vector<std::pair<DoseTypeGy, VoxelGridID> >::iterator pairItMax = maximumPositions->begin();
 				count = 0;
 
-				for (; pairIt2 != spMyResultPairs2->end() && count < 100; ++pairIt2) //output max. 100 maximum
+				for (; pairItMax != maximumPositions->end() && count < 100; ++pairItMax) //output max. 100 maximum
 				{
-					VoxelGridID voxelID = pairIt2->second;
+					VoxelGridID voxelID = pairItMax->second;
 					pt.add(statisticsTag + ".maximum.voxelGridID", voxelID);
 
-					VoxelGridIndex3D voxelIndex3D;
-					spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D);
-					std::string voxelIndex3DStr;
-					voxelIndex3DStr = boost::lexical_cast<std::string>(voxelIndex3D.x()) + ","
-					                  + boost::lexical_cast<std::string>(voxelIndex3D.y()) + ","
-					                  + boost::lexical_cast<std::string>(voxelIndex3D.z());
-					pt.add(statisticsTag + ".maximum.voxelIndex3D", voxelIndex3DStr);
+					//VoxelGridIndex3D voxelIndex3D;
+					//spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D);
+					//std::string voxelIndex3DStr;
+					//voxelIndex3DStr = boost::lexical_cast<std::string>(voxelIndex3D.x()) + ","
+					//                  + boost::lexical_cast<std::string>(voxelIndex3D.y()) + ","
+					//                  + boost::lexical_cast<std::string>(voxelIndex3D.z());
+					//pt.add(statisticsTag + ".maximum.voxelIndex3D", voxelIndex3DStr);
 
 
-					WorldCoordinate3D worldCoor;
-					spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor);
-					std::string worldCoorStr;
-					worldCoorStr = boost::lexical_cast<std::string>(worldCoor.x()) + ","
-					               + boost::lexical_cast<std::string>(worldCoor.y()) + ","
-					               + boost::lexical_cast<std::string>(worldCoor.z());
-					pt.add(statisticsTag + ".maximum.worldCoordinate", worldCoorStr);
+					//WorldCoordinate3D worldCoor;
+					//spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor);
+					//std::string worldCoorStr;
+					//worldCoorStr = boost::lexical_cast<std::string>(worldCoor.x()) + ","
+					//               + boost::lexical_cast<std::string>(worldCoor.y()) + ","
+					//               + boost::lexical_cast<std::string>(worldCoor.z());
+					//pt.add(statisticsTag + ".maximum.worldCoordinate", worldCoorStr);
 
 					count ++;
 				}
 
 				pt.put(statisticsTag + ".mean", aDoseStatistics->getMean());
 				pt.put(statisticsTag + ".standardDeviation", aDoseStatistics->getStdDeviation());
 				pt.put(statisticsTag + ".variance", aDoseStatistics->getVariance());
 
-				/*to do: x should be defined based on the user's feedback*/
-				//Dx
-				boost::property_tree::ptree dxNode1;
-				boost::property_tree::ptree dxNode2;
-				boost::property_tree::ptree dxNode3;
-				boost::property_tree::ptree dxNode4;
-				boost::property_tree::ptree dxNode5;
-				boost::property_tree::ptree dxNode6;
-
-				double absoluteVolume = aDoseStatistics->getVx(0);
-
-				dxNode1.put("", aDoseStatistics->getDx(absoluteVolume * 0.02));
-				dxNode2.put("", aDoseStatistics->getDx(absoluteVolume * 0.05));
-				dxNode3.put("", aDoseStatistics->getDx(absoluteVolume * 0.10));
-				dxNode4.put("", aDoseStatistics->getDx(absoluteVolume * 0.90));
-				dxNode5.put("", aDoseStatistics->getDx(absoluteVolume * 0.95));
-				dxNode6.put("", aDoseStatistics->getDx(absoluteVolume * 0.98));
-
-				dxNode1.put(xmlattrTag, 2);
-				dxNode2.put(xmlattrTag, 5);
-				dxNode3.put(xmlattrTag, 10);
-				dxNode4.put(xmlattrTag, 90);
-				dxNode5.put(xmlattrTag, 95);
-				dxNode6.put(xmlattrTag, 98);
-
-				pt.add_child(statisticsTag + ".Dx", dxNode1);
-				pt.add_child(statisticsTag + ".Dx", dxNode2);
-				pt.add_child(statisticsTag + ".Dx", dxNode3);
-				pt.add_child(statisticsTag + ".Dx", dxNode4);
-				pt.add_child(statisticsTag + ".Dx", dxNode5);
-				pt.add_child(statisticsTag + ".Dx", dxNode6);
+				double absoluteVolume = aDoseStatistics->getVolume();
+
+				try
+				{
+					//Dx
+					boost::property_tree::ptree dxNode1;
+					boost::property_tree::ptree dxNode2;
+					boost::property_tree::ptree dxNode3;
+					boost::property_tree::ptree dxNode4;
+					boost::property_tree::ptree dxNode5;
+					boost::property_tree::ptree dxNode6;
+
+					dxNode1.put("", aDoseStatistics->getDx(absoluteVolume * 0.02));
+					dxNode2.put("", aDoseStatistics->getDx(absoluteVolume * 0.05));
+					dxNode3.put("", aDoseStatistics->getDx(absoluteVolume * 0.10));
+					dxNode4.put("", aDoseStatistics->getDx(absoluteVolume * 0.90));
+					dxNode5.put("", aDoseStatistics->getDx(absoluteVolume * 0.95));
+					dxNode6.put("", aDoseStatistics->getDx(absoluteVolume * 0.98));
+
+					dxNode1.put(xmlattrTag, 2);
+					dxNode2.put(xmlattrTag, 5);
+					dxNode3.put(xmlattrTag, 10);
+					dxNode4.put(xmlattrTag, 90);
+					dxNode5.put(xmlattrTag, 95);
+					dxNode6.put(xmlattrTag, 98);
+
+					pt.add_child(statisticsTag + ".Dx", dxNode1);
+					pt.add_child(statisticsTag + ".Dx", dxNode2);
+					pt.add_child(statisticsTag + ".Dx", dxNode3);
+					pt.add_child(statisticsTag + ".Dx", dxNode4);
+					pt.add_child(statisticsTag + ".Dx", dxNode5);
+					pt.add_child(statisticsTag + ".Dx", dxNode6);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
 
 				if (aReferenceDose <= 0)
 				{
 					throw core::InvalidParameterException("aReferenceDose should be >0!");
 				}
 
-				//Vx
-				boost::property_tree::ptree vxNode1;
-				boost::property_tree::ptree vxNode2;
-				boost::property_tree::ptree vxNode3;
-				boost::property_tree::ptree vxNode4;
-				boost::property_tree::ptree vxNode5;
-				boost::property_tree::ptree vxNode6;
-
-				vxNode1.put("", aDoseStatistics->getVx(aReferenceDose * 0.02));
-				vxNode2.put("", aDoseStatistics->getVx(aReferenceDose * 0.05));
-				vxNode3.put("", aDoseStatistics->getVx(aReferenceDose * 0.10));
-				vxNode4.put("", aDoseStatistics->getVx(aReferenceDose * 0.90));
-				vxNode5.put("", aDoseStatistics->getVx(aReferenceDose * 0.95));
-				vxNode6.put("", aDoseStatistics->getVx(aReferenceDose * 0.98));
-
-				vxNode1.put(xmlattrTag, 2);
-				vxNode2.put(xmlattrTag, 5);
-				vxNode3.put(xmlattrTag, 10);
-				vxNode4.put(xmlattrTag, 90);
-				vxNode5.put(xmlattrTag, 95);
-				vxNode6.put(xmlattrTag, 98);
-
-				pt.add_child(statisticsTag + ".Vx", vxNode1);
-				pt.add_child(statisticsTag + ".Vx", vxNode2);
-				pt.add_child(statisticsTag + ".Vx", vxNode3);
-				pt.add_child(statisticsTag + ".Vx", vxNode4);
-				pt.add_child(statisticsTag + ".Vx", vxNode5);
-				pt.add_child(statisticsTag + ".Vx", vxNode6);
-
-				//MOHx
-				boost::property_tree::ptree mohxNode1;
-				boost::property_tree::ptree mohxNode2;
-				boost::property_tree::ptree mohxNode3;
-
-				mohxNode1.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.02));
-				mohxNode2.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.05));
-				mohxNode3.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.10));
-
-				mohxNode1.put(xmlattrTag, 2);
-				mohxNode2.put(xmlattrTag, 5);
-				mohxNode3.put(xmlattrTag, 10);
-
-				pt.add_child(statisticsTag + ".MOHx", mohxNode1);
-				pt.add_child(statisticsTag + ".MOHx", mohxNode2);
-				pt.add_child(statisticsTag + ".MOHx", mohxNode3);
-
-				//MOCx
-				boost::property_tree::ptree mocxNode1;
-				boost::property_tree::ptree mocxNode2;
-				boost::property_tree::ptree mocxNode3;
-
-				mocxNode1.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.02));
-				mocxNode2.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.05));
-				mocxNode3.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.10));
-
-				mocxNode1.put(xmlattrTag, 2);
-				mocxNode2.put(xmlattrTag, 5);
-				mocxNode3.put(xmlattrTag, 10);
-
-				pt.add_child(statisticsTag + ".MOCx", mocxNode1);
-				pt.add_child(statisticsTag + ".MOCx", mocxNode2);
-				pt.add_child(statisticsTag + ".MOCx", mocxNode3);
-
-
-				//MaxOHx
-				boost::property_tree::ptree maxohxNode1;
-				boost::property_tree::ptree maxohxNode2;
-				boost::property_tree::ptree maxohxNode3;
-
-				maxohxNode1.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.02));
-				maxohxNode2.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.05));
-				maxohxNode3.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.10));
-
-				maxohxNode1.put(xmlattrTag, 2);
-				maxohxNode2.put(xmlattrTag, 5);
-				maxohxNode3.put(xmlattrTag, 10);
-
-				pt.add_child(statisticsTag + ".MaxOHx", maxohxNode1);
-				pt.add_child(statisticsTag + ".MaxOHx", maxohxNode2);
-				pt.add_child(statisticsTag + ".MaxOHx", maxohxNode3);
-
-				//MinOCx
-				boost::property_tree::ptree minocxNode1;
-				boost::property_tree::ptree minocxNode2;
-				boost::property_tree::ptree minocxNode3;
-
-				minocxNode1.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.02));
-				minocxNode2.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.05));
-				minocxNode3.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.10));
-
-				minocxNode1.put(xmlattrTag, 2);
-				minocxNode2.put(xmlattrTag, 5);
-				minocxNode3.put(xmlattrTag, 10);
-
-				pt.add_child(statisticsTag + ".MinOCx", minocxNode1);
-				pt.add_child(statisticsTag + ".MinOCx", minocxNode2);
-				pt.add_child(statisticsTag + ".MinOCx", minocxNode3);
+				try
+				{
+					//Vx
+					boost::property_tree::ptree vxNode1;
+					boost::property_tree::ptree vxNode2;
+					boost::property_tree::ptree vxNode3;
+					boost::property_tree::ptree vxNode4;
+					boost::property_tree::ptree vxNode5;
+					boost::property_tree::ptree vxNode6;
+
+					vxNode1.put("", aDoseStatistics->getVx(aReferenceDose * 0.02));
+					vxNode2.put("", aDoseStatistics->getVx(aReferenceDose * 0.05));
+					vxNode3.put("", aDoseStatistics->getVx(aReferenceDose * 0.10));
+					vxNode4.put("", aDoseStatistics->getVx(aReferenceDose * 0.90));
+					vxNode5.put("", aDoseStatistics->getVx(aReferenceDose * 0.95));
+					vxNode6.put("", aDoseStatistics->getVx(aReferenceDose * 0.98));
+
+					vxNode1.put(xmlattrTag, 2);
+					vxNode2.put(xmlattrTag, 5);
+					vxNode3.put(xmlattrTag, 10);
+					vxNode4.put(xmlattrTag, 90);
+					vxNode5.put(xmlattrTag, 95);
+					vxNode6.put(xmlattrTag, 98);
+
+					pt.add_child(statisticsTag + ".Vx", vxNode1);
+					pt.add_child(statisticsTag + ".Vx", vxNode2);
+					pt.add_child(statisticsTag + ".Vx", vxNode3);
+					pt.add_child(statisticsTag + ".Vx", vxNode4);
+					pt.add_child(statisticsTag + ".Vx", vxNode5);
+					pt.add_child(statisticsTag + ".Vx", vxNode6);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
+
+				try
+				{
+					//MOHx
+					boost::property_tree::ptree mohxNode1;
+					boost::property_tree::ptree mohxNode2;
+					boost::property_tree::ptree mohxNode3;
+
+					mohxNode1.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.02));
+					mohxNode2.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.05));
+					mohxNode3.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.10));
+
+					mohxNode1.put(xmlattrTag, 2);
+					mohxNode2.put(xmlattrTag, 5);
+					mohxNode3.put(xmlattrTag, 10);
+
+					pt.add_child(statisticsTag + ".MOHx", mohxNode1);
+					pt.add_child(statisticsTag + ".MOHx", mohxNode2);
+					pt.add_child(statisticsTag + ".MOHx", mohxNode3);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
+
+				try
+				{
+					//MOCx
+					boost::property_tree::ptree mocxNode1;
+					boost::property_tree::ptree mocxNode2;
+					boost::property_tree::ptree mocxNode3;
+
+					mocxNode1.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.02));
+					mocxNode2.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.05));
+					mocxNode3.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.10));
+
+					mocxNode1.put(xmlattrTag, 2);
+					mocxNode2.put(xmlattrTag, 5);
+					mocxNode3.put(xmlattrTag, 10);
+
+					pt.add_child(statisticsTag + ".MOCx", mocxNode1);
+					pt.add_child(statisticsTag + ".MOCx", mocxNode2);
+					pt.add_child(statisticsTag + ".MOCx", mocxNode3);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
+
+				try
+				{
+					//MaxOHx
+					boost::property_tree::ptree maxohxNode1;
+					boost::property_tree::ptree maxohxNode2;
+					boost::property_tree::ptree maxohxNode3;
+
+					maxohxNode1.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.02));
+					maxohxNode2.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.05));
+					maxohxNode3.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.10));
+
+					maxohxNode1.put(xmlattrTag, 2);
+					maxohxNode2.put(xmlattrTag, 5);
+					maxohxNode3.put(xmlattrTag, 10);
+
+					pt.add_child(statisticsTag + ".MaxOHx", maxohxNode1);
+					pt.add_child(statisticsTag + ".MaxOHx", maxohxNode2);
+					pt.add_child(statisticsTag + ".MaxOHx", maxohxNode3);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
+
+				try
+				{
+					//MinOCx
+					boost::property_tree::ptree minocxNode1;
+					boost::property_tree::ptree minocxNode2;
+					boost::property_tree::ptree minocxNode3;
+
+					minocxNode1.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.02));
+					minocxNode2.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.05));
+					minocxNode3.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.10));
+
+					minocxNode1.put(xmlattrTag, 2);
+					minocxNode2.put(xmlattrTag, 5);
+					minocxNode3.put(xmlattrTag, 10);
+
+					pt.add_child(statisticsTag + ".MinOCx", minocxNode1);
+					pt.add_child(statisticsTag + ".MinOCx", minocxNode2);
+					pt.add_child(statisticsTag + ".MinOCx", minocxNode3);
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
 				return pt;
 
 			}
 
 			void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName,
 			                         DoseTypeGy aReferenceDose)
 			{
 				boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose);
 
 				try
 				{
 					boost::property_tree::xml_parser::write_xml(aFileName, pt);
 				}
 				catch (boost::property_tree::xml_parser_error& /*e*/)
 				{
 					throw core::InvalidParameterException("Write xml failed: xml_parser_error!");
 				}
 			}
 
 			XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose)
 			{
 				boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose);
 				std::stringstream sstr;
 
 				try
 				{
 					boost::property_tree::xml_parser::write_xml(sstr, pt);
 				}
 				catch (boost::property_tree::xml_parser_error& /*e*/)
 				{
 					throw core::InvalidParameterException("Write xml to string failed: xml_parser_error!");
 				}
 
 				return sstr.str();
 			}
 
 
 			StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics,
 			        DoseTypeGy aReferenceDose)
 			{
 
 				std::stringstream sstr;
 
 
-				rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator =
-				    aDoseStatistics->getDoseIterator();
-				rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor =
-				    spDoseIterator->getDoseAccessor();
+				//rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator =
+				//    aDoseStatistics->getDoseIterator();
+				//rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor =
+				//    spDoseIterator->getDoseAccessor();
 
-				boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs =
-				    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
-				ResultListPointer spMyResultPairs(myResultPairs);
-				boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs2 =
-				    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
-				ResultListPointer spMyResultPairs2(myResultPairs2);
+				//boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs =
+				//    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+				//ResultListPointer spMyResultPairs(myResultPairs);
+				//boost::shared_ptr< std::vector<std::pair<DoseTypeGy, VoxelGridID> > > myResultPairs2 =
+				//    boost::make_shared< std::vector<std::pair<DoseTypeGy, VoxelGridID> > >();
+				//ResultListPointer spMyResultPairs2(myResultPairs2);
 
-				sstr << aDoseStatistics->getVx(0) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getMaximum(spMyResultPairs2) << columnSeparator;
-				sstr << aDoseStatistics->getMinimum(spMyResultPairs) << columnSeparator;
+				sstr << aDoseStatistics->getVolume() * 1000 << columnSeparator; // cm3 to mm3
+				sstr << aDoseStatistics->getMaximum() << columnSeparator;
+				sstr << aDoseStatistics->getMinimum() << columnSeparator;
 
 				sstr << aDoseStatistics->getMean() << columnSeparator;
 				sstr << aDoseStatistics->getStdDeviation() << columnSeparator;
 				sstr << aDoseStatistics->getVariance() << columnSeparator;
 
 				/*to do: x should be defined based on the user's feedback*/
 				if (aReferenceDose <= 0)
 				{
 					throw core::InvalidParameterException("aReferenceDose should be >0!");
 				}
 
 				//Vx
 
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.02) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.05) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.10) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.90) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.95) * 1000 << columnSeparator; // cm3 to mm3
-				sstr << aDoseStatistics->getVx(aReferenceDose * 0.98) * 1000 << columnSeparator; // cm3 to mm3
+				try
+				{
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.02) * 1000 << columnSeparator; // cm3 to mm3
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.05) * 1000 << columnSeparator; // cm3 to mm3
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.10) * 1000 << columnSeparator; // cm3 to mm3
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.90) * 1000 << columnSeparator; // cm3 to mm3
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.95) * 1000 << columnSeparator; // cm3 to mm3
+					sstr << aDoseStatistics->getVx(aReferenceDose * 0.98) * 1000 << columnSeparator; // cm3 to mm3
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
 				//Dx
 
-				double absoluteVolume = aDoseStatistics->getVx(0);
-
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.02) << columnSeparator;
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.05) << columnSeparator;
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.10) << columnSeparator;
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.90) << columnSeparator;
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.95) << columnSeparator;
-				sstr << aDoseStatistics->getDx(absoluteVolume * 0.98) << columnSeparator;
-
-
-				//MOHx
-				sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.02) << columnSeparator;
-				sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.05) << columnSeparator;
-				sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.10) << columnSeparator;
+				double absoluteVolume = aDoseStatistics->getVolume();
 
+				try
+				{
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.02) << columnSeparator;
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.05) << columnSeparator;
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.10) << columnSeparator;
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.90) << columnSeparator;
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.95) << columnSeparator;
+					sstr << aDoseStatistics->getDx(absoluteVolume * 0.98) << columnSeparator;
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
-				//MOCx
-				sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.02) << columnSeparator;
-				sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.05) << columnSeparator;
-				sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.10) << columnSeparator;
+				try
+				{
+					//MOHx
+					sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.02) << columnSeparator;
+					sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.05) << columnSeparator;
+					sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.10) << columnSeparator;
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
+				try
+				{
+					//MOCx
+					sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.02) << columnSeparator;
+					sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.05) << columnSeparator;
+					sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.10) << columnSeparator;
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
-				//MaxOHx
-				sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.02) << columnSeparator;
-				sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.05) << columnSeparator;
-				sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.10) << columnSeparator;
+				try
+				{
+					//MaxOHx
+					sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.02) << columnSeparator;
+					sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.05) << columnSeparator;
+					sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.10) << columnSeparator;
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
-				//MinOCx
-				sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.02) << columnSeparator;
-				sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.05) << columnSeparator;
-				sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.10) << columnSeparator;
+				try
+				{
+					//MinOCx
+					sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.02) << columnSeparator;
+					sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.05) << columnSeparator;
+					sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.10) << columnSeparator;
+				}
+				catch (core::DataNotAvailableException)
+				{
+					//as data is not available (was not computed by doseStatistics), it cannot be written
+				}
 
 
 				return sstr.str();
 
 			}
 
 		}
 	}
 }
\ No newline at end of file
diff --git a/testing/algorithms/CMakeLists.txt b/testing/algorithms/CMakeLists.txt
index 4a106e5..6d4bc45 100644
--- a/testing/algorithms/CMakeLists.txt
+++ b/testing/algorithms/CMakeLists.txt
@@ -1,21 +1,22 @@
 #-----------------------------------------------------------------------------
 # Setup the system information test.  Write out some basic failsafe
 # information in case the test doesn't run.
 #-----------------------------------------------------------------------------
 
 
 SET(ALGORITHMS_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbAlgorithmsTests)
 SET(ALGORITHMS_HEADER_TEST ${EXECUTABLE_OUTPUT_PATH}/rttbAlgorithmsHeaderTest)
 
 SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data)
 
 SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary)
 
 
 #-----------------------------------------------------------------------------
 ADD_TEST(DoseStatisticsTest ${ALGORITHMS_TESTS} DoseStatisticsTest)
+ADD_TEST(DoseStatisticsCalculatorTest ${ALGORITHMS_TESTS} DoseStatisticsCalculatorTest)
 ADD_TEST(ArithmeticTest ${ALGORITHMS_TESTS} ArithmeticTest)
 ADD_TEST(BinaryFunctorDoseAccessorTest ${ALGORITHMS_TESTS} BinaryFunctorDoseAccessorTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/dicompylerTestDose.dcm")
 
 RTTB_CREATE_TEST_MODULE(rttbAlgorithms DEPENDS RTTBAlgorithms RTTBMasks RTTBDicomIO PACKAGE_DEPENDS Boost Litmus DCMTK)
 
diff --git a/testing/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp
new file mode 100644
index 0000000..507ce7c
--- /dev/null
+++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp
@@ -0,0 +1,314 @@
+// -----------------------------------------------------------------------
+// RTToolbox - DKFZ radiotherapy quantitative evaluation library
+//
+// Copyright (c) German Cancer Research Center (DKFZ),
+// Software development for Integrated Diagnostics and Therapy (SIDT).
+// ALL RIGHTS RESERVED.
+// See rttbCopyright.txt or
+// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
+//
+// This software is distributed WITHOUT ANY WARRANTY; without even
+// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+// PURPOSE.  See the above copyright notices for more information.
+//
+//------------------------------------------------------------------------
+/*!
+// @file
+// @version $Revision: 707 $ (last changed revision)
+// @date    $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date)
+// @author  $Author: floca $ (last changed by)
+*/
+
+#include <boost/make_shared.hpp>
+#include <boost/shared_ptr.hpp>
+
+#include "litCheckMacros.h"
+
+#include "rttbBaseType.h"
+#include "rttbGenericDoseIterator.h"
+#include "rttbDoseIteratorInterface.h"
+#include "rttbNullPointerException.h"
+#include "rttbDoseStatisticsCalculator.h"
+#include "rttbInvalidDoseException.h"
+#include "rttbInvalidParameterException.h"
+#include "rttbDataNotAvailableException.h"
+
+#include "../core/DummyDoseAccessor.h"
+
+namespace rttb
+{
+	namespace testing
+	{
+
+		typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer;
+		typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer;
+		typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer;
+		typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer;
+
+		/*! @brief DoseStatisticsCalculatorTest - test the API of DoseStatisticsCalculator
+			1) test constructors
+			2) test setDoseIterator
+			3) test calculateDoseSatistics
+			4) get statistical values
+		*/
+
+		int DoseStatisticsCalculatorTest(int argc, char* argv[])
+		{
+			PREPARE_DEFAULT_TEST_REPORTING;
+
+			boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor = boost::make_shared<DummyDoseAccessor>();
+			DoseAccessorPointer spDoseAccessor(spTestDoseAccessor);
+			const std::vector<DoseTypeGy>* doseVals = spTestDoseAccessor->getDoseVector();
+
+			boost::shared_ptr<core::GenericDoseIterator> spTestDoseIterator =
+			    boost::make_shared<core::GenericDoseIterator>(spDoseAccessor);
+			DoseIteratorPointer spDoseIterator(spTestDoseIterator);
+
+			DoseIteratorPointer spDoseIteratorNull;
+
+			//1) test constructors
+			// the values cannot be accessed from outside, therefore correct default values are not tested
+			CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur);
+			rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur;
+			rttb::algorithms::DoseStatisticsCalculator myDoseStatCalculaturToBeFilled;
+
+			CHECK_THROW_EXPLICIT(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIteratorNull),
+			                     core::NullPointerException);
+
+			CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator));
+			rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator);
+
+			//2) test setDoseIterator
+			CHECK_THROW_EXPLICIT(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIteratorNull), core::NullPointerException);
+			CHECK_NO_THROW(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIterator));
+
+			//3) test calculateDoseStatistics
+			DoseStatisticsPointer theStatistics, theStatistics2;
+			CHECK_THROW_EXPLICIT(myEmptyDoseStatCalculatur.calculateDoseStatistics(), core::NullPointerException);
+			CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics());
+			CHECK_NO_THROW(theStatistics2 = myDoseStatCalculaturToBeFilled.calculateDoseStatistics());
+			CHECK_EQUAL(theStatistics->getMinimumPositions()->empty(), true);
+			CHECK_EQUAL(theStatistics->getMaximumPositions()->empty(), true);
+			CHECK_EQUAL(theStatistics->getAllVx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllDx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllVx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllMaxOHx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllMOHx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllMOCx().empty(), true);
+			CHECK_EQUAL(theStatistics->getAllMinOCx().empty(), true);
+
+			//check default values for computeComplexMeasures=true
+			DoseStatisticsPointer theStatisticsDefault;
+			CHECK_NO_THROW(theStatisticsDefault = myDoseStatsCalculator.calculateDoseStatistics(true));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.02 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.05 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.1 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.2 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.8 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.9 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.95 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getVx(0.98 * theStatisticsDefault->getMaximum()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.02 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.05 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.1 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.2 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.8 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.9 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.95 * theStatisticsDefault->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatisticsDefault->getDx(0.98 * theStatisticsDefault->getNumberOfVoxels()));
+
+			//check manually set values for omputeComplexMeasures=true
+			std::vector<double> precomputeDoseValues, precomputeVolumeValues;
+			precomputeDoseValues.push_back(0.01 * theStatistics->getMaximum());
+			precomputeDoseValues.push_back(0.02 * theStatistics->getMaximum());
+			precomputeDoseValues.push_back(0.05 * theStatistics->getMaximum());
+
+			precomputeVolumeValues.push_back(0.9 * theStatistics->getNumberOfVoxels());
+			precomputeVolumeValues.push_back(0.95 * theStatistics->getNumberOfVoxels());
+			precomputeVolumeValues.push_back(0.99 * theStatistics->getNumberOfVoxels());
+
+			CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(true, precomputeDoseValues,
+			                               precomputeVolumeValues));
+			CHECK_NO_THROW(theStatistics->getVx(0.01 * theStatistics->getMaximum()));
+			CHECK_NO_THROW(theStatistics->getVx(0.02 * theStatistics->getMaximum()));
+			CHECK_NO_THROW(theStatistics->getVx(0.05 * theStatistics->getMaximum()));
+			CHECK_THROW_EXPLICIT(theStatistics->getVx(0.03 * theStatistics->getMaximum()), core::DataNotAvailableException);
+			CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels()));
+			CHECK_NO_THROW(theStatistics->getDx(0.99 * theStatistics->getNumberOfVoxels()));
+			CHECK_THROW_EXPLICIT(theStatistics->getDx(0.03 * theStatistics->getNumberOfVoxels()), core::DataNotAvailableException);
+			CHECK_EQUAL(theStatistics->getVx(0.02 * theStatistics->getMaximum()),
+			            theStatisticsDefault->getVx(0.02 * theStatistics->getMaximum()));
+			CHECK_EQUAL(theStatistics->getVx(0.05 * theStatistics->getMaximum()),
+			            theStatisticsDefault->getVx(0.05 * theStatistics->getMaximum()));
+			CHECK_EQUAL(theStatistics->getDx(0.9 * theStatistics->getNumberOfVoxels()),
+			            theStatisticsDefault->getDx(0.9 * theStatistics->getNumberOfVoxels()));
+			CHECK_EQUAL(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels()),
+			            theStatisticsDefault->getDx(0.95 * theStatistics->getNumberOfVoxels()));
+
+			//MOHx, MOCx, MaxOHx and MinOCx are computed analogous to Dx, they will not be checked.
+
+			//4) get statistical values
+			CHECK_EQUAL(theStatistics->getNumberOfVoxels(), doseVals->size());
+			CHECK_EQUAL(theStatistics2->getNumberOfVoxels(), doseVals->size());
+
+			//compute simple statistical values (min, mean, max, stddev) for comparison
+			DoseStatisticType maximum = 0;
+			DoseStatisticType minimum = 1000000;
+			DoseStatisticType mean = 0;
+			DoseStatisticType variance = 0;
+			std::vector<DoseTypeGy>::const_iterator doseIt = doseVals->begin();
+
+			while (doseIt != doseVals->end())
+			{
+				if (maximum < *doseIt)
+				{
+					maximum = *doseIt;
+				}
+
+				if (minimum > *doseIt)
+				{
+					minimum = *doseIt;
+				}
+
+				mean += *doseIt;
+				++doseIt;
+			}
+
+			mean /= doseVals->size();
+			float compMean = (int(mean * 100)) / 100;
+
+			doseIt = doseVals->begin();
+
+			while (doseIt != doseVals->end())
+			{
+				variance += pow(*doseIt - mean, 2);
+				++doseIt;
+			}
+
+			variance /= doseVals->size() - 1;
+			DoseStatisticType stdDev = pow(variance, 0.5);
+
+			CHECK_EQUAL(theStatistics->getMaximum(), maximum);
+			CHECK_EQUAL(theStatistics2->getMaximum(), maximum);
+
+			CHECK_EQUAL(theStatistics->getMinimum(), minimum);
+			CHECK_EQUAL(theStatistics2->getMinimum(), minimum);
+
+			CHECK_EQUAL(theStatistics->getMean(), mean);
+			CHECK_EQUAL(theStatistics2->getMean(), mean);
+
+			CHECK_EQUAL(theStatistics->getStdDeviation(), stdDev);
+			CHECK_EQUAL(theStatistics2->getStdDeviation(), stdDev);
+
+			CHECK_EQUAL(theStatistics->getVariance(), variance);
+			CHECK_EQUAL(theStatistics2->getVariance(), variance);
+
+			//check for complex doseStatistics (maximumPositions, minimumPositions, Vx, Dx, MOHx, MOCx, MAXOHx, MinOCx)
+
+			unsigned int nMax = 0, nMin = 0;
+
+			doseIt = doseVals->begin();
+
+			while (doseIt != doseVals->end())
+			{
+				if (*doseIt == theStatistics->getMaximum())
+				{
+					nMax++;
+				}
+
+				if (*doseIt == theStatistics->getMinimum())
+				{
+					nMin++;
+				}
+
+				++doseIt;
+			}
+
+			//only 100 positions are stored
+			if (nMax > 100)
+			{
+				nMax = 100;
+			}
+
+			if (nMin > 100)
+			{
+				nMin = 100;
+			}
+
+			auto maximaPositions = theStatistics->getMaximumPositions();
+			auto minimaPositions = theStatistics->getMinimumPositions();
+
+			CHECK_EQUAL(maximaPositions->size(), nMax);
+			CHECK_EQUAL(minimaPositions->size(), nMin);
+
+			for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions);
+			     ++maximaPositionsIterator)
+			{
+				CHECK_EQUAL(maximaPositionsIterator->first, theStatistics->getMaximum());
+			}
+
+			for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions);
+			     ++minimaPositionsIterator)
+			{
+				CHECK_EQUAL(minimaPositionsIterator->first, theStatistics->getMinimum());
+			}
+
+
+			//generate specific example dose
+			maximum = 9.5;
+			minimum = 2.5;
+			mean = 6;
+			int sizeTemplate = 500;
+			std::vector<DoseTypeGy> aDoseVector;
+
+			for (int i = 0; i < sizeTemplate; i++)
+			{
+				aDoseVector.push_back(maximum);
+				aDoseVector.push_back(minimum);
+			}
+
+			core::GeometricInfo geoInfo = spTestDoseAccessor->getGeometricInfo();
+			geoInfo.setNumRows(20);
+			geoInfo.setNumColumns(10);
+			geoInfo.setNumSlices(5);
+
+			boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor2 =
+			    boost::make_shared<DummyDoseAccessor>(aDoseVector, geoInfo);
+			DoseAccessorPointer spDoseAccessor2(spTestDoseAccessor2);
+
+			boost::shared_ptr<core::GenericDoseIterator> spTestDoseIterator2 =
+			    boost::make_shared<core::GenericDoseIterator>(spDoseAccessor2);
+			DoseIteratorPointer spDoseIterator2(spTestDoseIterator2);
+
+			rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator2(spDoseIterator2);
+			DoseStatisticsPointer theStatistics3 = myDoseStatsCalculator2.calculateDoseStatistics();
+
+			CHECK_EQUAL(theStatistics3->getMaximum(), maximum);
+			CHECK_EQUAL(theStatistics3->getMinimum(), minimum);
+			CHECK_EQUAL(theStatistics3->getMean(), mean);
+
+
+			maximaPositions = theStatistics3->getMaximumPositions();
+			minimaPositions = theStatistics3->getMinimumPositions();
+
+			CHECK_EQUAL(maximaPositions->empty(), false);
+			CHECK_EQUAL(minimaPositions->empty(), false);
+
+			for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions);
+			     ++maximaPositionsIterator)
+			{
+				CHECK_EQUAL(maximaPositionsIterator->first, theStatistics3->getMaximum());
+			}
+
+			for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions);
+			     ++minimaPositionsIterator)
+			{
+				CHECK_EQUAL(minimaPositionsIterator->first, theStatistics3->getMinimum());
+			}
+
+			RETURN_AND_REPORT_TEST_SUCCESS;
+		}
+
+	}//end namespace testing
+}//end namespace rttb
diff --git a/testing/algorithms/DoseStatisticsTest.cpp b/testing/algorithms/DoseStatisticsTest.cpp
index 5c6a803..10bc265 100644
--- a/testing/algorithms/DoseStatisticsTest.cpp
+++ b/testing/algorithms/DoseStatisticsTest.cpp
@@ -1,217 +1,218 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE.  See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date    $Date$ (last change date)
 // @author  $Author$ (last changed by)
 */
 
 #include <boost/make_shared.hpp>
 #include <boost/shared_ptr.hpp>
 
 #include "litCheckMacros.h"
 
 #include "rttbBaseType.h"
-#include "rttbGenericDoseIterator.h"
-#include "rttbDoseIteratorInterface.h"
-#include "rttbNullPointerException.h"
 #include "rttbDoseStatistics.h"
-#include "rttbInvalidDoseException.h"
+#include "rttbDataNotAvailableException.h"
 
-#include "../core/DummyDoseAccessor.h"
+namespace rttb
+{
+	namespace testing
+	{
 
-namespace rttb{
-	namespace testing{
-
-		typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer;
-		typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer;
 		typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer;
+		typedef rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType DoseToVolumeFunctionType;
+		typedef rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType;
 
 		/*! @brief DoseStatisticsTest - test the API of DoseStatistics
-			1) test constructors 
-			2) test setDoseIterator
-			3) test getNumberOfVoxels
-			4) get statistical values
-			floating point accuracy requires this otherwise the test 4 fails!
+			1) test constructors
+			2) test setters
+			3) test getters of complex statistics (with stored key and without stored key)
 		*/
 
-		int DoseStatisticsTest(int argc, char* argv[] )
-			{
+		int DoseStatisticsTest(int argc, char* argv[])
+		{
 			PREPARE_DEFAULT_TEST_REPORTING;
 
-			boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor = boost::make_shared<DummyDoseAccessor>();
-			DoseAccessorPointer spDoseAccessor(spTestDoseAccessor);
-			const std::vector<DoseTypeGy>* doseVals = spTestDoseAccessor->getDoseVector();
-
-			boost::shared_ptr<core::GenericDoseIterator> spTestDoseIterator = 
-				boost::make_shared<core::GenericDoseIterator>(spDoseAccessor);
-			DoseIteratorPointer spDoseIterator (spTestDoseIterator);
-
-			DoseIteratorPointer spDoseIteratorNull;
-
-			//1) test constructors 
-			// the values cannot be accessed from outside, therefore correct default values are not tested
-			CHECK_NO_THROW(rttb::algorithms::DoseStatistics myEmptyDoseStat);
-			rttb::algorithms::DoseStatistics myEmptyDoseStat;
-			rttb::algorithms::DoseStatistics myDoseStatToBeFilled;
-
-			CHECK_THROW_EXPLICIT(rttb::algorithms::DoseStatistics myDoseStats(spDoseIteratorNull),core::NullPointerException);
-
-			CHECK_NO_THROW(rttb::algorithms::DoseStatistics myDoseStats(spDoseIterator));
-			rttb::algorithms::DoseStatistics myDoseStats(spDoseIterator);
-
-			//2) test setDoseIterator
-			CHECK_THROW_EXPLICIT(myDoseStatToBeFilled.setDoseIterator(spDoseIteratorNull),core::NullPointerException);
-			CHECK_NO_THROW(myDoseStatToBeFilled.setDoseIterator(spDoseIterator));
-
-			//3) test getNumberOfVoxels
-			CHECK_THROW_EXPLICIT(myEmptyDoseStat.getNumberOfVoxels(),core::InvalidDoseException);
-			CHECK_EQUAL(myDoseStats.getNumberOfVoxels(),doseVals->size());
-			CHECK_EQUAL(myDoseStatToBeFilled.getNumberOfVoxels(),doseVals->size());	  
-
-			//4) get statistical values
-      CHECK_THROW_EXPLICIT(myEmptyDoseStat.getNumberOfVoxels(),core::InvalidDoseException);
-			CHECK_EQUAL(myDoseStats.getNumberOfVoxels(),doseVals->size());
-			CHECK_EQUAL(myDoseStatToBeFilled.getNumberOfVoxels(),doseVals->size());	
-
-			//compute values for comparison
-			DoseStatisticType maximum = 0;
-			DoseStatisticType minimum = 1000000;
-			DoseStatisticType mean = 0;
-			std::vector<DoseTypeGy>::const_iterator doseIt = doseVals->begin(); 
-			while(doseIt != doseVals->end())
-				{
-				if (maximum < *doseIt)
-					{
-					maximum = *doseIt;
-					}
-				if (minimum > *doseIt)
-					{
-					minimum = *doseIt;
-					}
-				mean += *doseIt;
-				doseIt++;
-				}
-			mean /=doseVals->size();
-			float compMean = (int(mean*100))/100;
-
-			boost::shared_ptr< std::vector<std::pair<DoseTypeGy,VoxelGridID> > > myResultPairs = 
-				boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-			ResultListPointer spMyResultPairs(myResultPairs);
-			boost::shared_ptr< std::vector<std::pair<DoseTypeGy,VoxelGridID> > > myResultPairs2 = 
-				boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-			ResultListPointer spMyResultPairs2(myResultPairs2);
-
-      CHECK_THROW_EXPLICIT(myEmptyDoseStat.getMaximum(spMyResultPairs),core::InvalidDoseException);
-			CHECK_EQUAL(myDoseStats.getMaximum(spMyResultPairs),maximum);
-			CHECK_EQUAL(myDoseStatToBeFilled.getMaximum(spMyResultPairs2),maximum);
-			CHECK(spMyResultPairs->size()>0);
-
-			CHECK_EQUAL(spMyResultPairs->size(),spMyResultPairs2->size());
-			std::vector<std::pair<DoseTypeGy,VoxelGridID> >::iterator pairIt = spMyResultPairs->begin();
-			std::vector<std::pair<DoseTypeGy,VoxelGridID> >::iterator pairIt2 = spMyResultPairs2->begin();
-
-			for (;pairIt != spMyResultPairs->end();++pairIt)
-				{
-				CHECK_EQUAL(pairIt->first,pairIt2->first);
-				CHECK_EQUAL(pairIt->second,pairIt2->second);
-				++pairIt2;
-				}
-
-			spMyResultPairs = boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-			spMyResultPairs2 = boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-      CHECK_THROW_EXPLICIT(myEmptyDoseStat.getMinimum(spMyResultPairs),core::InvalidDoseException);
-			CHECK_EQUAL(myDoseStats.getMinimum(spMyResultPairs),minimum);
-			CHECK_EQUAL(myDoseStatToBeFilled.getMinimum(spMyResultPairs2),minimum);
-			CHECK(spMyResultPairs->size()>0);
-
-			CHECK_EQUAL(spMyResultPairs->size(),spMyResultPairs2->size());
-			pairIt = spMyResultPairs->begin();
-			pairIt2 = spMyResultPairs2->begin();
-
-			for (;pairIt != spMyResultPairs->end();++pairIt)
-				{
-				CHECK_EQUAL(pairIt->first,pairIt2->first);
-				CHECK_EQUAL(pairIt->second,pairIt2->second);
-				++pairIt2;
-				}
-
-      CHECK_THROW_EXPLICIT(myEmptyDoseStat.getMean(),core::InvalidDoseException);
-			float tmpMean = myDoseStats.getMean();
-			tmpMean = (int(tmpMean*100))/100;
-			CHECK_EQUAL(tmpMean,compMean);
-			tmpMean = myDoseStatToBeFilled.getMean();
-			tmpMean = (int(tmpMean*100))/100;
-			CHECK_EQUAL(tmpMean,compMean);
-
-			//generate specific example dose
-			maximum = 9.5;
-			minimum = 2.5;
-			mean = 6;
-			int sizeTemplate = 500;
-			std::vector<DoseTypeGy> aDoseVector;
-			for (int i = 0; i < sizeTemplate; i++){	
-				aDoseVector.push_back(maximum);
-				aDoseVector.push_back(minimum);
-				}
-
-			core::GeometricInfo geoInfo = spTestDoseAccessor->getGeometricInfo();
-			geoInfo.setNumRows(20);
-			geoInfo.setNumColumns(10);
-			geoInfo.setNumSlices(5);
-
-			boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor2 = 
-				boost::make_shared<DummyDoseAccessor>(aDoseVector,geoInfo);
-			DoseAccessorPointer spDoseAccessor2(spTestDoseAccessor2);
-
-			boost::shared_ptr<core::GenericDoseIterator> spTestDoseIterator2 = 
-				boost::make_shared<core::GenericDoseIterator>(spDoseAccessor2);
-			DoseIteratorPointer spDoseIterator2 (spTestDoseIterator2);
-
-			rttb::algorithms::DoseStatistics myDoseStats2(spDoseIterator2);
-
-			spMyResultPairs = boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-			CHECK_EQUAL(myDoseStats2.getMaximum(spMyResultPairs),maximum);
-			CHECK(spMyResultPairs->size()>0);
-
-			CHECK_EQUAL(spMyResultPairs->size(),sizeTemplate);
-			pairIt = spMyResultPairs->begin();
-
-			for (;pairIt != spMyResultPairs->end();++pairIt)
-				{
-				CHECK_EQUAL(pairIt->first,maximum);
-				}
-
-			spMyResultPairs = boost::make_shared< std::vector<std::pair<DoseTypeGy,VoxelGridID> > >();
-			int maxmialFound = 200; 
-			CHECK_EQUAL(myDoseStats2.getMinimum(spMyResultPairs, maxmialFound),minimum);
-			CHECK(spMyResultPairs->size()>0);
-
-			CHECK_EQUAL(spMyResultPairs->size(),maxmialFound);
-			pairIt = spMyResultPairs->begin();
-
-			for (;pairIt != spMyResultPairs->end();++pairIt)
-				{
-				CHECK_EQUAL(pairIt->first,minimum);
-				}
-
-			CHECK_EQUAL(myDoseStats2.getMean(),mean);
+			DoseStatisticType minimum = 1.0;
+			DoseStatisticType mean = 5.5;
+			DoseStatisticType maximum = 108.2;
+			DoseStatisticType stdDeviation = 10.1;
+			unsigned int numVoxels = 100000;
+			VolumeType volume = numVoxels * (0.5 * 0.5 * 0.5);
+
+			std::vector<std::pair<DoseTypeGy, VoxelGridID> > minVoxels;
+			std::vector<std::pair<DoseTypeGy, VoxelGridID> > maxVoxels;
+
+			minVoxels.push_back(std::make_pair(1.0, 11));
+			minVoxels.push_back(std::make_pair(1.0, 22));
+			minVoxels.push_back(std::make_pair(1.0, 33));
+			minVoxels.push_back(std::make_pair(1.0, 44));
+
+			maxVoxels.push_back(std::make_pair(108.2, 5));
+			maxVoxels.push_back(std::make_pair(108.2, 6));
+			maxVoxels.push_back(std::make_pair(108.2, 7));
+			maxVoxels.push_back(std::make_pair(108.2, 8));
+
+			ResultListPointer resultsMinVoxels = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >(minVoxels);
+			ResultListPointer resultsMaxVoxels = boost::make_shared<std::vector<std::pair<DoseTypeGy, VoxelGridID> > >(maxVoxels);
+
+			DoseToVolumeFunctionType Vx;
+			Vx.insert(std::make_pair(1.1, 1000));
+			Vx.insert(std::make_pair(106.9, 99000));
+
+			VolumeToDoseFunctionType Dx;
+			Dx.insert(std::make_pair(1000, 1.1));
+			Dx.insert(std::make_pair(99000, 106.9));
+
+			VolumeToDoseFunctionType MOHx;
+			MOHx.insert(std::make_pair(1000, 5));
+			MOHx.insert(std::make_pair(99000, 105.5));
+
+			VolumeToDoseFunctionType MOCx;
+			MOCx.insert(std::make_pair(1000, 10));
+			MOCx.insert(std::make_pair(99000, 99));
+
+			VolumeToDoseFunctionType MaxOHx;
+			MaxOHx.insert(std::make_pair(1000, 40));
+			MaxOHx.insert(std::make_pair(99000, 98.3));
+
+			VolumeToDoseFunctionType MinOCx;
+			MinOCx.insert(std::make_pair(1000, 25.5));
+			MinOCx.insert(std::make_pair(99000, 102.7));
+
+
+			//1) test constructors
+
+			CHECK_NO_THROW(rttb::algorithms::DoseStatistics aDoseStatistic(minimum, maximum, mean, stdDeviation, numVoxels,
+			               volume));
+
+			rttb::algorithms::DoseStatistics aDoseStatistic(minimum, maximum, mean, stdDeviation, numVoxels, volume);
+			CHECK_EQUAL(aDoseStatistic.getMinimum(), minimum);
+			CHECK_EQUAL(aDoseStatistic.getMaximum(), maximum);
+			CHECK_EQUAL(aDoseStatistic.getMean(), mean);
+			CHECK_EQUAL(aDoseStatistic.getStdDeviation(), stdDeviation);
+			CHECK_EQUAL(aDoseStatistic.getVariance(), stdDeviation * stdDeviation);
+			CHECK_EQUAL(aDoseStatistic.getNumberOfVoxels(), numVoxels);
+			CHECK_EQUAL(aDoseStatistic.getVolume(), volume);
+
+			//check default values for unset complex values
+			CHECK_EQUAL(aDoseStatistic.getMaximumPositions()->empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getMinimumPositions()->empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllDx().empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllVx().empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllMOHx().empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllMOCx().empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllMaxOHx().empty(), true);
+			CHECK_EQUAL(aDoseStatistic.getAllMinOCx().empty(), true);
+
+			CHECK_NO_THROW(rttb::algorithms::DoseStatistics aDoseStatisticComplex(minimum, maximum, mean, stdDeviation, numVoxels,
+			               volume,
+			               resultsMaxVoxels, resultsMinVoxels, Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx));
+			rttb::algorithms::DoseStatistics aDoseStatisticComplex(minimum, maximum, mean, stdDeviation, numVoxels, volume,
+			        resultsMaxVoxels, resultsMinVoxels, Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx);
+
+			CHECK_EQUAL(aDoseStatisticComplex.getMaximumPositions(), resultsMaxVoxels);
+			CHECK_EQUAL(aDoseStatisticComplex.getMinimumPositions(), resultsMinVoxels);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllDx() == Dx, true);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllVx() == Vx, true);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllMOHx() == MOHx, true);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllMOCx() == MOCx, true);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllMaxOHx() == MaxOHx, true);
+			CHECK_EQUAL(aDoseStatisticComplex.getAllMinOCx() == MinOCx, true);
+
+			//2) test setters (only complex statistics have setters)
+			CHECK_NO_THROW(aDoseStatistic.setMaximumVoxelPositions(resultsMaxVoxels));
+			CHECK_NO_THROW(aDoseStatistic.setMinimumVoxelPositions(resultsMinVoxels));
+			CHECK_NO_THROW(aDoseStatistic.setDx(Dx));
+			CHECK_NO_THROW(aDoseStatistic.setVx(Vx));
+			CHECK_NO_THROW(aDoseStatistic.setMOHx(MOHx));
+			CHECK_NO_THROW(aDoseStatistic.setMOCx(MOCx));
+			CHECK_NO_THROW(aDoseStatistic.setMaxOHx(MaxOHx));
+			CHECK_NO_THROW(aDoseStatistic.setMinOCx(MinOCx));
+
+			CHECK_EQUAL(aDoseStatistic.getMaximumPositions(), resultsMaxVoxels);
+			CHECK_EQUAL(aDoseStatistic.getMinimumPositions(), resultsMinVoxels);
+			CHECK_EQUAL(aDoseStatistic.getAllDx() == Dx, true);
+			CHECK_EQUAL(aDoseStatistic.getAllVx() == Vx, true);
+			CHECK_EQUAL(aDoseStatistic.getAllMOHx() == MOHx, true);
+			CHECK_EQUAL(aDoseStatistic.getAllMOCx() == MOCx, true);
+			CHECK_EQUAL(aDoseStatistic.getAllMaxOHx() == MaxOHx, true);
+			CHECK_EQUAL(aDoseStatistic.getAllMinOCx() == MinOCx, true);
+
+			//3) test getters of complex statistics(with stored key and without stored key)
+			//getAll*() already tested in (2)
+			Vx.clear();
+			Vx.insert(std::make_pair(1.1, 1000));
+			Vx.insert(std::make_pair(5.0, 2300));
+			Vx.insert(std::make_pair(90, 90500));
+			Vx.insert(std::make_pair(107, 99000));
+
+			Dx.clear();
+			Dx.insert(std::make_pair(1000, 1.1));
+			Dx.insert(std::make_pair(2000, 2.0));
+			Dx.insert(std::make_pair(5000, 10.8));
+			Dx.insert(std::make_pair(90000, 89.5));
+			Dx.insert(std::make_pair(98000, 104.4));
+			Dx.insert(std::make_pair(99000, 106.9));
+
+			rttb::algorithms::DoseStatistics aDoseStatisticNewValues(minimum, maximum, mean, stdDeviation, numVoxels, volume);
+			aDoseStatisticNewValues.setDx(Dx);
+			aDoseStatisticNewValues.setVx(Vx);
+
+			CHECK_NO_THROW(aDoseStatisticNewValues.getVx(1.1));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getVx(90));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getDx(1000));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getDx(98000));
+
+			CHECK_EQUAL(aDoseStatisticNewValues.getVx(1.1), Vx.find(1.1)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getVx(90), Vx.find(90)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getDx(1000), Dx.find(1000)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getDx(98000), Dx.find(98000)->second);
+
+			//test if key-value combination NOT in map
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getDx(1001), core::DataNotAvailableException);
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getVx(10), core::DataNotAvailableException);
+
+			double closestDxKey, closestVxKey;
+			CHECK_NO_THROW(aDoseStatisticNewValues.getDx(900, true, closestDxKey));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getDx(99001, true, closestDxKey));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getVx(10, true, closestVxKey));
+			CHECK_EQUAL(aDoseStatisticNewValues.getDx(900, true, closestDxKey), Dx.find(1000)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getDx(99001, true, closestDxKey), Dx.find(99000)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getVx(10, true, closestVxKey), Vx.find(5.0)->second);
+			CHECK_EQUAL(closestDxKey, 99000);
+			CHECK_EQUAL(closestVxKey, 5);
+
+			//equal distance to two values. First value is returned.
+			CHECK_NO_THROW(aDoseStatisticNewValues.getDx(1500, true, closestDxKey));
+			CHECK_NO_THROW(aDoseStatisticNewValues.getVx(98.5, true, closestVxKey));
+			CHECK_EQUAL(aDoseStatisticNewValues.getDx(1500, true, closestDxKey), Dx.find(1000)->second);
+			CHECK_EQUAL(aDoseStatisticNewValues.getVx(98.5, true, closestVxKey), Vx.find(90.0)->second);
+			CHECK_EQUAL(closestDxKey, 1000);
+			CHECK_EQUAL(closestVxKey, 90.0);
+
+			double dummy;
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx(25), core::DataNotAvailableException);
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx(9999), core::DataNotAvailableException);
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx(25, true, dummy), core::DataNotAvailableException);
+			CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx(9999, true, dummy), core::DataNotAvailableException);
 
 			RETURN_AND_REPORT_TEST_SUCCESS;
-			}
+		}
 
-		}//end namespace testing
-	}//end namespace rttb
+	}//end namespace testing
+}//end namespace rttb
diff --git a/testing/algorithms/files.cmake b/testing/algorithms/files.cmake
index 31ddd28..e056ae8 100644
--- a/testing/algorithms/files.cmake
+++ b/testing/algorithms/files.cmake
@@ -1,19 +1,20 @@
 SET(CPP_FILES 
 	DoseStatisticsTest.cpp
+	DoseStatisticsCalculatorTest.cpp
 	ArithmeticTest.cpp
 	BinaryFunctorDoseAccessorTest.cpp
 	rttbAlgorithmsTests.cpp
 	
 	#include dummy accessor files
 	../core/DummyDoseAccessor.cpp
 	../core/DummyMaskAccessor.cpp
 	../core/DummyMutableDoseAccessor.cpp
 	
   )
 
 SET(H_FILES 
 	#include dummy accessor files
 	../core/DummyDoseAccessor.h
 	../core/DummyMaskAccessor.h
 	../core/DummyMutableDoseAccessor.h
 )
diff --git a/testing/algorithms/rttbAlgorithmsTests.cpp b/testing/algorithms/rttbAlgorithmsTests.cpp
index e03c3e9..69263f6 100644
--- a/testing/algorithms/rttbAlgorithmsTests.cpp
+++ b/testing/algorithms/rttbAlgorithmsTests.cpp
@@ -1,65 +1,66 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE.  See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date    $Date$ (last change date)
 // @author  $Author$ (last changed by)
 */
 
 // this file defines the rttbAlgorithmsTests for the test driver
 // and all it expects is that you have a function called RegisterTests
 #if defined(_MSC_VER)
 #pragma warning ( disable : 4786 )
 #endif
 
 
 #include "litMultiTestsMain.h"
 
 namespace rttb
 {
 	namespace testing
 	{
 
 		void registerTests()
 		{
 			LIT_REGISTER_TEST(DoseStatisticsTest);
+			LIT_REGISTER_TEST(DoseStatisticsCalculatorTest);
 			LIT_REGISTER_TEST(ArithmeticTest);
 			LIT_REGISTER_TEST(BinaryFunctorDoseAccessorTest);
 		}
 	}
 }
 
 int main(int argc, char* argv[])
 {
 	int result = 0;
 
 	rttb::testing::registerTests();
 
 	try
 	{
 		result = lit::multiTestsMain(argc, argv);
 	}
 	catch (const std::exception& /*e*/)
 	{
 		result = -1;
 	}
 	catch (...)
 	{
 		result = -1;
 	}
 
 	return result;
 }
diff --git a/testing/io/other/DoseStatisticsIOTest.cpp b/testing/io/other/DoseStatisticsIOTest.cpp
index 4fe32e1..c76beb0 100644
--- a/testing/io/other/DoseStatisticsIOTest.cpp
+++ b/testing/io/other/DoseStatisticsIOTest.cpp
@@ -1,98 +1,106 @@
 // -----------------------------------------------------------------------
 // RTToolbox - DKFZ radiotherapy quantitative evaluation library
 //
 // Copyright (c) German Cancer Research Center (DKFZ),
 // Software development for Integrated Diagnostics and Therapy (SIDT).
 // ALL RIGHTS RESERVED.
 // See rttbCopyright.txt or
 // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^]
 //
 // This software is distributed WITHOUT ANY WARRANTY; without even
 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 // PURPOSE. See the above copyright notices for more information.
 //
 //------------------------------------------------------------------------
 /*!
 // @file
 // @version $Revision$ (last changed revision)
 // @date $Date$ (last change date)
 // @author $Author$ (last changed by)
 */
 
 // this file defines the rttbCoreTests for the test driver
 // and all it expects is that you have a function called RegisterTests
 
 #include <boost/make_shared.hpp>
 #include <boost/shared_ptr.hpp>
 #include <boost/lexical_cast.hpp>
 
 #include "litCheckMacros.h"
 
 #include "rttbBaseType.h"
-#include "rttbDoseStatistics.h"
+#include "rttbDoseStatisticsCalculator.h"
 #include "rttbDoseStatisticsXMLWriter.h"
 #include "rttbGenericDoseIterator.h"
 #include "rttbDoseIteratorInterface.h"
 #include "rttbInvalidParameterException.h"
 #include "rttbNullPointerException.h"
 
 #include "../../core/DummyDoseAccessor.h"
 
 
 namespace rttb
 {
 
 	namespace testing
 	{
 
 		/*! @brief OtherIOTest - test the IO for dose statistics
 		    1) test exception
-			2) test writing statistcs to xml file
+			2) test writing statistics to xml file
 		*/
 
 		int DoseStatisticsIOTest(int argc, char* argv[])
 		{
 			typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer;
 			typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer;
-			typedef boost::shared_ptr<rttb::algorithms::DoseStatistics> DoseStatisticsPtr;
+			typedef boost::shared_ptr<rttb::algorithms::DoseStatisticsCalculator> DoseStatisticsCalculatorPtr;
 
 			PREPARE_DEFAULT_TEST_REPORTING;
 
 			/* generate dummy dose */
 			boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor = boost::make_shared<DummyDoseAccessor>();
 			DoseAccessorPointer spDoseAccessor(spTestDoseAccessor);
 
 			boost::shared_ptr<core::GenericDoseIterator> spTestDoseIterator =
 			    boost::make_shared<core::GenericDoseIterator>(spDoseAccessor);
 			DoseIteratorPointer spDoseIterator(spTestDoseIterator);
 
-			rttb::algorithms::DoseStatistics myDoseStats(spDoseIterator);
-			DoseStatisticsPtr myDoseStatsPtr = boost::make_shared<rttb::algorithms::DoseStatistics>
-			                                   (myDoseStats);
+			rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator);
+			auto myDoseStatsSimple = myDoseStatsCalculator.calculateDoseStatistics();
+			auto myDoseStatsComplex = myDoseStatsCalculator.calculateDoseStatistics(true);
 
 			/* test exception */
-			CHECK_THROW_EXPLICIT(io::other::writeDoseStatistics(myDoseStatsPtr, "test.test", 0),
+			CHECK_THROW_EXPLICIT(io::other::writeDoseStatistics(myDoseStatsSimple, "test.test", 0),
 			                     core::InvalidParameterException);
 
 
-			/* test writing statistcs to xml file */
-			FileNameString fN = "testStatistics.xml";
-			CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsPtr, fN));
+			/* test writing statistics to xml file */
+			FileNameString filenameSimple = "testStatisticsSimple.xml";
+			CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsSimple, filenameSimple));
 
-			/* test writing statistcs to string */
-			boost::property_tree::ptree pt = io::other::writeDoseStatistics(myDoseStatsPtr);
-			XMLString str = io::other::writerDoseStatisticsToString(myDoseStatsPtr);
+			FileNameString filenameComplex = "testStatisticsComplex.xml";
+			CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsComplex, filenameComplex));
 
-			std::stringstream sstr;
-			boost::property_tree::xml_parser::write_xml(sstr, pt);
-			CHECK_EQUAL(str, sstr.str());
+			/* test writing statistics to string */
+			boost::property_tree::ptree ptSimple = io::other::writeDoseStatistics(myDoseStatsSimple);
+			XMLString strSimple = io::other::writerDoseStatisticsToString(myDoseStatsSimple);
 
+			boost::property_tree::ptree ptComplex = io::other::writeDoseStatistics(myDoseStatsComplex);
+			XMLString strComplex = io::other::writerDoseStatisticsToString(myDoseStatsComplex);
 
+			std::stringstream sstrSimple;
+			boost::property_tree::xml_parser::write_xml(sstrSimple, ptSimple);
+			CHECK_EQUAL(strSimple, sstrSimple.str());
+
+			std::stringstream sstrComplex;
+			boost::property_tree::xml_parser::write_xml(sstrComplex, ptComplex);
+			CHECK_EQUAL(strComplex, sstrComplex.str());
 
 
 			RETURN_AND_REPORT_TEST_SUCCESS;
 		}
 
 	}//testing
 }//rttb