diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp
index 390b772..10a1bf2 100644
--- a/code/masks/boost/rttbBoostMask.cpp
+++ b/code/masks/boost/rttbBoostMask.cpp
@@ -1,269 +1,269 @@
 #include <iostream>
 #include <deque>
 #include <algorithm>
 
 #include <boost/geometry/geometries/register/point.hpp>
 #include <boost/geometry/geometries/register/ring.hpp>
 #include <boost/geometry/io/wkt/wkt.hpp>
 #include <boost/geometry/multi/geometries/multi_polygon.hpp>
 #include <boost/make_shared.hpp>
 
 #include "rttbBoostMask.h"
 #include "rttbNullPointerException.h"
 #include "rttbInvalidParameterException.h"
 
 namespace rttb
 {
 	namespace masks
 	{
 		namespace boost
 		{
 			BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure)
 				:_geometricInfo(aDoseGeoInfo), _structure(aStructure), _voxelInStructure(::boost::make_shared<MaskVoxelList>()){
 
 					_isUpToDate = false;
 					if(! _geometricInfo ){
 						throw rttb::core::NullPointerException("Error: Geometric info is NULL!");
 					}
 					else if(!_structure){
 						throw rttb::core::NullPointerException("Error: Structure is NULL!");
 					}
 
 			}
 
 			void BoostMask::calcMask(){
 				if(hasSelfIntersections()){
 					throw rttb::core::InvalidParameterException("Error: structure has self intersections!");
 				}
 
 				for(unsigned int i=0; i< _geometricInfo->getNumSlices(); i++){
-					VoxelIndexVector voxelList = getBoundingBox(i);
+					BoostMask::BoostRingVector intersectionSlicePolygons;
+					VoxelIndexVector voxelList = getBoundingBox(i, intersectionSlicePolygons);
 					rttb::VoxelGridIndex3D gridIndex3D0 = voxelList.at(0);
 					rttb::VoxelGridIndex3D gridIndex3D1 = voxelList.at(1);
 
 
 					for(unsigned int indexX = gridIndex3D0[0] ; indexX <= gridIndex3D1[0]; indexX ++  ){
 						for(unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY ++){
 							rttb::VoxelGridIndex3D currentIndex;
 							currentIndex[0] = indexX;
 							currentIndex[1] = indexY;
 							currentIndex[2] = i;
 
 							rttb::VoxelGridID gridID; 
 							_geometricInfo->convert(currentIndex, gridID);
 
-							std::deque<BoostMask::BoostPolygon2D> polygons = getIntersections(currentIndex);
+							std::deque<BoostMask::BoostPolygon2D> polygons = getIntersections(currentIndex, intersectionSlicePolygons);
 							double intersectionArea = calcArea(polygons);
 
 							double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1];
 
 							if(intersectionArea > 0) {
 								double volumeFraction = intersectionArea/voxelSize;
 								if(volumeFraction > 1 && (volumeFraction-1) <= 0.0001){
 									volumeFraction = 1;
 								}
 								core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, (volumeFraction));
 								_voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure
 							}
 						}
 					}
 				}
 
 				sort(_voxelInStructure->begin(), _voxelInStructure->end());
 			}
 
 			bool BoostMask::hasSelfIntersections(){
 				bool hasSelfIntersection = false;
 
 				for(unsigned int indexZ=0; indexZ< _geometricInfo->getNumSlices(); indexZ++){
 					rttb::VoxelGridIndex3D currentIndex(0,0,indexZ);
 
-					BoostMask::BoostRingVector boostPolygons = getIntersectionSlicePolygons(currentIndex);
+					BoostMask::BoostRingVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]);
 
 					for(unsigned int i =0; i<boostPolygons.size(); i++){
 						::boost::geometry::correct(boostPolygons.at(i));
 						//test if polygon has self intersection
 						if(::boost::geometry::detail::overlay::has_self_intersections(boostPolygons.at(i))){
 							hasSelfIntersection = true;
 							std::cout << "Has self intersection! Slice: "<< indexZ << ", polygon "<<i <<std::endl;
 							break;
 						}
 						//test if the polygons on the same slice has intersection
 						for(unsigned int j =0; j<boostPolygons.size() && j!=i; j++){
 							::boost::geometry::correct(boostPolygons.at(j));
 							BoostPolygonDeque intersection;
 							::boost::geometry::intersection(boostPolygons.at(i), boostPolygons.at(j), intersection);
 							if(intersection.size()>0){
 								hasSelfIntersection = true;
 								std::cout << "Has self intersection! Slice: "<< indexZ << ", polygon "<<i << " and "<<j<<std::endl;
 								break;
 							}
 						}
 					}	
 
 				}
 
 				return hasSelfIntersection;
 			}
 
 
 			/*Get the 4 voxel index of the bounding box of the polygon in the slice with sliceNumber*/
-			BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber){
+			BoostMask::VoxelIndexVector BoostMask::getBoundingBox(const unsigned int sliceNumber, BoostRingVector& intersectionSlicePolygons){
 				if(sliceNumber<0 || sliceNumber >= _geometricInfo->getNumSlices()){
 					throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!"));
 				}
 
 				rttb::VoxelGridIndex3D currentIndex(0,0,sliceNumber);
 
-				BoostRingVector boostPolygons = getIntersectionSlicePolygons(currentIndex);
+				intersectionSlicePolygons = getIntersectionSlicePolygons(currentIndex[2]);
 
 				double xMin = 0;
 				double yMin = 0;
 				double xMax = 0; 
 				double yMax = 0;
-				for(unsigned int i =0; i<boostPolygons.size(); i++){
-					::boost::geometry::correct(boostPolygons.at(i));
+				for(unsigned int i =0; i<intersectionSlicePolygons.size(); i++){
+					::boost::geometry::correct(intersectionSlicePolygons.at(i));
 					::boost::geometry::model::box<BoostPoint2D> box;
-					::boost::geometry::envelope(boostPolygons.at(i), box);
+					::boost::geometry::envelope(intersectionSlicePolygons.at(i), box);
 					BoostPoint2D minPoint = box.min_corner();
 					BoostPoint2D maxPoint = box.max_corner();
 					if(i==0){
 						xMin = minPoint.x();
 						yMin = minPoint.y();
 						xMax = maxPoint.x();
 						yMax = maxPoint.y();
 					}
 					xMin = std::min(xMin, minPoint.x());
 					yMin = std::min(yMin, minPoint.y());
 					xMax = std::max(xMax, maxPoint.x());
 					yMax = std::max(yMax, maxPoint.y());
 				}
 
 				rttb::WorldCoordinate3D nullWorldCoord;
 				_geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0,0,0),nullWorldCoord);
 				rttb::WorldCoordinate3D minWorldCoord(xMin, yMin, nullWorldCoord.z());
 				rttb::WorldCoordinate3D maxWorldCoord(xMax, yMax, nullWorldCoord.z());
 
 				rttb::VoxelGridIndex3D index0;
 				rttb::VoxelGridIndex3D index1;
 				_geometricInfo->worldCoordinateToIndex(minWorldCoord, index0);
 				_geometricInfo->worldCoordinateToIndex(maxWorldCoord, index1);
 
 				VoxelIndexVector voxelList;
 				voxelList.push_back(index0);
 				voxelList.push_back(index1);
 
 				return voxelList;
 			} 
 
 
 			/*Get intersection polygons of the contour and a voxel polygon*/
-			BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D){
+			BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostRingVector& intersectionSlicePolygons){
 				BoostMask::BoostPolygonDeque polygonDeque;
 
 				BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D);
 				::boost::geometry::correct(voxelPolygon);
 
-				BoostMask::BoostRingVector contourVector = getIntersectionSlicePolygons(aVoxelIndex3D);
-				for(unsigned int i=0; i<contourVector.size(); i++){
-					BoostRing2D contour = contourVector.at(i);
+				for(unsigned int i=0; i<intersectionSlicePolygons.size(); i++){
+					BoostRing2D contour = intersectionSlicePolygons.at(i);
 					::boost::geometry::correct(contour);
 					BoostPolygonDeque intersection;
 					::boost::geometry::intersection(voxelPolygon, contour, intersection);
 
 					polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end());
 				}
 
 				return polygonDeque;
 			}
 
 
 			/*Calculate the intersection area*/
 			double BoostMask::calcArea(std::deque<BoostMask::BoostPolygon2D> aPolygonDeque){
 				double area = 0;
 				for(unsigned int i=0; i<aPolygonDeque.size(); i++){
 					area += ::boost::geometry::area(aPolygonDeque.at(i));
 				}
 				return area;
 			}
 
 			VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel){
 				rttb::VoxelGridIndex3D gridIndex3D; 
 				_geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D);
 				return gridIndex3D;
 			}
 
 			BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector(){
 				if(!_isUpToDate){
 					calcMask();
 				}
 				return _voxelInStructure;
 			}
 
 			BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon){
 				BoostMask::BoostRing2D polygon2D;
 				BoostPoint2D firstPoint;
 				for(unsigned int i=0; i<aRTTBPolygon.size(); i++){
 					rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i);
 					BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); 
 					if(i==0){
 						firstPoint = boostPoint;
 					}
 					::boost::geometry::append(polygon2D, boostPoint);
 				}
 				::boost::geometry::append(polygon2D, firstPoint);
 				return polygon2D;
 
 			}
 
-			BoostMask::BoostRingVector BoostMask::getIntersectionSlicePolygons(const rttb::VoxelGridIndex3D& aVoxelGrid3D){
+			BoostMask::BoostRingVector BoostMask::getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ){
 				BoostRingVector boostPolygonVector;
 
 				rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector();
 
 				for(unsigned int i=0; i<polygonSequence.size(); i++){
 					PolygonType rttbPolygon = polygonSequence.at(i);
 					if(rttbPolygon.size()>0){
 						double polygonZCoor = rttbPolygon.at(0)[2];
 						rttb::WorldCoordinate3D polygonPoint(0,0, polygonZCoor);
 						rttb::VoxelGridIndex3D polygonPointIndex3D;
 						_geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D);
 
 						//if the z
-						if(aVoxelGrid3D[2] == polygonPointIndex3D[2]){
+						if(aVoxelGridIndexZ == polygonPointIndex3D[2]){
 							boostPolygonVector.push_back(convert(rttbPolygon));
 						}
 					}
 				}
 
 				return boostPolygonVector;
 			}
 
 			BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D){
 				BoostMask::BoostRing2D polygon;
 				rttb::WorldCoordinate3D rttbPoint;
 				_geometricInfo->indexToWorldCoordinate(aVoxelGrid3D, rttbPoint);
 
 				BoostPoint2D point1 (rttbPoint[0], rttbPoint[1]);
 				::boost::geometry::append(polygon, point1);
 
 				double xSpacing = _geometricInfo->getSpacing()[0];
 				double ySpacing = _geometricInfo->getSpacing()[1];
 
 				BoostPoint2D point2(rttbPoint[0]+xSpacing, rttbPoint[1]);
 				::boost::geometry::append(polygon, point2);
 
 				BoostPoint2D point3(rttbPoint[0]+xSpacing, rttbPoint[1]+ySpacing);
 				::boost::geometry::append(polygon, point3);
 
 				BoostPoint2D point4(rttbPoint[0], rttbPoint[1]+ySpacing);
 				::boost::geometry::append(polygon, point4);
 
 				::boost::geometry::append(polygon, point1);
 
 				return polygon;
 
 			}
 		}
 	}
 }
\ No newline at end of file
diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/boost/rttbBoostMask.h
index d908ed7..9e8f16c 100644
--- a/code/masks/boost/rttbBoostMask.h
+++ b/code/masks/boost/rttbBoostMask.h
@@ -1,118 +1,118 @@
 // -----------------------------------------------------------------------
 // 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 __BOOST_MASK_H
 #define __BOOST_MASK_H
 
 #include "rttbBaseType.h"
 #include "rttbStructure.h"
 #include "rttbGeometricInfo.h"
 #include "rttbMaskVoxel.h"
 #include "rttbMaskAccessorInterface.h"
 
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/point_xy.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
 #include <boost/geometry/geometries/ring.hpp>
 #include <boost/shared_ptr.hpp>
 
 namespace rttb
 {
 	namespace masks
 	{
 		namespace boost
 		{
 			class BoostMask
 			{	
 
 			public:
 				typedef ::boost::shared_ptr<rttb::core::GeometricInfo> GeometricInfoPointer;
 				typedef core::Structure::StructTypePointer StructPointer;
 				typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList;
 				typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer;
 
 				/*! @brief Constructor
 				* @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL
 				*/
 				BoostMask(GeometricInfoPointer aDoseGeoInfo, StructPointer aStructure);
 
 				/*! @brief Generate mask and return the voxels in the mask
 				* @exception rttb::core::InvalidParameterException thrown if the structure has self intersections
 				*/
 				MaskVoxelListPointer getRelevantVoxelVector();
 
 			private:
 				typedef ::boost::geometry::model::d2::point_xy<double> BoostPoint2D;
 				typedef ::boost::geometry::model::polygon<::boost::geometry::model::d2::point_xy<double> > BoostPolygon2D;
 				typedef ::boost::geometry::model::ring<::boost::geometry::model::d2::point_xy<double> > BoostRing2D;
 				typedef std::deque<BoostPolygon2D> BoostPolygonDeque;
 				typedef std::vector<BoostRing2D> BoostRingVector;
 				typedef std::vector<rttb::VoxelGridIndex3D> VoxelIndexVector;
 
 				GeometricInfoPointer _geometricInfo;
 
 				StructPointer _structure;
 
 				//vector of the MaskVoxel inside the structure
 				MaskVoxelListPointer _voxelInStructure;
 
 				/*! @brief If the mask is up to date
 				*/
 				bool _isUpToDate;
 
 				/*! @brief Voxelization and generate mask
 				*/
 				void calcMask();
 
 				/*! @brief Check if the structure has self intersections*/
 				bool hasSelfIntersections();
 
 				/*! @brief Get the min/max voxel index of the bounding box of the polygon in the slice with sliceNumber*/
-				VoxelIndexVector getBoundingBox(unsigned int sliceNumber);
+				VoxelIndexVector getBoundingBox(const unsigned int sliceNumber, BoostRingVector& intersectionSlicePolygons);
 
 				/*! @brief Get intersection polygons of the contour and a voxel polygon*/
-				BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D);
+				BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostRingVector& intersectionSlicePolygons);
 
 				/*! @brief Calculate the intersection area*/
 				double calcArea(BoostPolygonDeque aPolygonDeque);
 
 				/*! @brief Get grid index of a mask voxel*/
 				VoxelGridIndex3D getGridIndex3D(const core::MaskVoxel& aMaskVoxel);
 
 				/*! @brief Convert RTTB polygon to boost polygon*/
 				BoostRing2D convert(const rttb::PolygonType& aRTTBPolygon);
 
 				/*! @brief Get the intersection slice of the voxel, return the polygons in the slice*/
-				BoostRingVector getIntersectionSlicePolygons(const rttb::VoxelGridIndex3D& aVoxelGrid3D);
+				BoostRingVector getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ);
 
 				/*! @brief Get the voxel 2d contour polygon*/
 				BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D);
 
 			};
 
 		}
 
 
 	}
 }
 
 #endif
\ No newline at end of file
diff --git a/testing/examples/CMakeLists.txt b/testing/examples/CMakeLists.txt
index fa2b62b..0ff90f3 100644
--- a/testing/examples/CMakeLists.txt
+++ b/testing/examples/CMakeLists.txt
@@ -1,40 +1,40 @@
 #-----------------------------------------------------------------------------
 # Setup the system information test.  Write out some basic failsafe
 # information in case the test doesn't run.
 #-----------------------------------------------------------------------------
 
 
 SET(CORE_TEST_EXAMPLES ${EXECUTABLE_OUTPUT_PATH}/rttbExamplesTests)
 
 SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data)
 
 SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary)
 
 
 #-----------------------------------------------------------------------------
 ADD_TEST(RTBioModelExampleTest ${CORE_TEST_EXAMPLES} RTBioModelExampleTest
 "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt"
 "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt" 
 "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk6.txt"
 "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk8.txt")
 ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest 
 "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling.dcm"
 "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling05.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantFiftyGridScaling01.dcm")
 ADD_TEST(RTDoseIndexTest ${CORE_TEST_EXAMPLES} RTDoseIndexTest 
 "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt"
 "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt")
 ADD_TEST(RTDoseStatisticsTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsTest 
 "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" 
 "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm")
 ADD_TEST(RTDVHTest ${CORE_TEST_EXAMPLES} RTDVHTest "${TEST_DATA_ROOT}/TestDVH/dvh_test.txt")
 ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest 
 "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm"
 "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm")
 ADD_TEST(RTBioModelScatterPlotExampleTest ${CORE_TEST_EXAMPLES} RTBioModelScatterPlotExampleTest
 "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt")
 ADD_TEST(VoxelizationValidationTest ${CORE_TEST_EXAMPLES} VoxelizationValidationTest 
-"${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" 
-"${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/boostMask/" "${TEST_DATA_ROOT}/OTBMask/")
+"${TEST_DATA_ROOT}/DICOM/StructureSet/H000068_20100316T175500_RTSTRUCT_676.IMA" 
+"${TEST_DATA_ROOT}/DICOM/TestDose/H000068_20100316T175500_RTDOSE_445.IMA" "${TEST_DATA_ROOT}/boostMask/" "${TEST_DATA_ROOT}/OTBMask/")
 
 RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBOtherIO RTTBITKIO RTTBModels PACKAGE_DEPENDS BoostBinaries Litmus DCMTK ITK)
 
diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp
index 6ef01fe..42358f2 100644
--- a/testing/examples/VoxelizationValidationTest.cpp
+++ b/testing/examples/VoxelizationValidationTest.cpp
@@ -1,161 +1,184 @@
 // -----------------------------------------------------------------------
 // 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: 929 $ (last changed revision)
 // @date $Date: 2015-04-08 14:50:57 +0200 (Mi, 08 Apr 2015) $ (last change date)
 // @author $Author: zhangl $ (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 <iomanip>
 
 #include "litCheckMacros.h"
 
 
 #include "dcmtk/config/osconfig.h"    /* make sure OS specific configuration is included first */
 
 #include "dcmtk/dcmrt/drtdose.h"
 #include "dcmtk/dcmrt/drtstrct.h"
 
 
 #include "rttbBaseType.h"
 #include "rttbDVHCalculator.h"
 #include "rttbGenericMaskedDoseIterator.h"
 #include "rttbGenericDoseIterator.h"
 #include "rttbNullPointerException.h"
 #include "rttbInvalidParameterException.h"
 #include "rttbDicomDoseAccessor.h"
 #include "rttbDicomFileDoseAccessorGenerator.h"
 #include "rttbDicomFileStructureSetGenerator.h"
 #include "rttbOTBMaskAccessor.h"
 #include "rttbDVHTxtFileReader.h"
 #include "rttbBoostMaskAccessor.h"
 #include "rttbITKImageMaskAccessorConverter.h"
 #include "rttbImageWriter.h"
 
 
 namespace rttb
 {
 
 	namespace testing
 	{
 
 		/*! @brief VoxelizationValidationTest.
 		
 		*/
 
 		int VoxelizationValidationTest(int argc, char* argv[])
 		{
 			typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer;
 			typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer;
 			typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer;
 			typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer;
 
 			PREPARE_DEFAULT_TEST_REPORTING;
 			//ARGUMENTS: 1: structure file name
 			//           2: dose1 file name
 			//           3: directory name to write boost mask of all structures
 			//           4: directory name to write OTB mask of all structures
 			
 
-			std::string RTSTRUCT_FILENAME;
-			std::string RTDOSE_FILENAME;
-			std::string BoostMask_DIRNAME;
-			std::string OTBMask_DIRNAME;
+			std::string RTSTRUCT_FILENAME = "D:/Temp/Test/str.IMA";
+			std::string RTDOSE_FILENAME= "D:/Temp/Test/dose.IMA";
+			std::string BoostMask_DIRNAME= "D:/Temp/Test/Boost/";
+			std::string OTBMask_DIRNAME= "D:/Temp/Test/OTB/";
 
 
-			if (argc > 4)
+			/*if (argc > 4)
 			{
 				RTSTRUCT_FILENAME = argv[1];
 				RTDOSE_FILENAME = argv[2];
 				BoostMask_DIRNAME = argv[3];
 				OTBMask_DIRNAME = argv[4];
-			}
+			}*/
 
 
 			OFCondition status;
 			DcmFileFormat fileformat;
 
 			/* read dicom-rt dose */
 			::DRTDoseIOD rtdose1;
 			io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str());
 			DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor());
 			boost::shared_ptr<core::GeometricInfo> geometricPtr = boost::make_shared<core::GeometricInfo>(doseAccessor1->getGeometricInfo());
 
 			//create a vector of MaskAccessors (one for each structure)
 			StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator(
 			        RTSTRUCT_FILENAME.c_str()).generateStructureSet();
 
 			std::vector<MaskAccessorPointer> rtStructSetMaskAccessorVec;
 			
 
 			if (rtStructureSet->getNumberOfStructures() > 0)
 			{
 				for (int j = 2; j < rtStructureSet->getNumberOfStructures(); j++)
 				{
 					clock_t start(clock());
 					//create OTB MaskAccessor
 					::boost::shared_ptr<masks::legacy::OTBMaskAccessor> spOTBMaskAccessor =
 					    ::boost::make_shared<masks::legacy::OTBMaskAccessor>(rtStructureSet->getStructure(j),
 					            doseAccessor1->getGeometricInfo());
 					spOTBMaskAccessor->updateMask();
 					MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor);
+					
+					::boost::shared_ptr<core::GenericMaskedDoseIterator> spMaskedDoseIteratorTmp =
+					    ::boost::make_shared<core::GenericMaskedDoseIterator>(spMaskAccessor, doseAccessor1);
+					DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp);
+					rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(),
+					                               doseAccessor1->getDoseUID());
+					rttb::core::DVH dvh = *(calc.generateDVH());
+
+					clock_t finish(clock());
+					std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl;
 					//write the mask image to a file	
 					rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor);
 					CHECK(itkConverter.process());
+					
 					std::stringstream fileNameSstr;
 					fileNameSstr<<OTBMask_DIRNAME<<j<<".mhd";
 					std::cout << fileNameSstr.str()<<std::endl;
 					rttb::io::itk::ImageWriter writer(fileNameSstr.str(), itkConverter.getITKImage());
 					CHECK(writer.writeITKFile());
-
-					clock_t finish(clock());
-					std::cout << "OTB Mask Calculation and write file time: " << finish - start << " ms" << std::endl;
+					
+					
 
 					clock_t start2(clock());
 					//create Boost MaskAccessor		
 					MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared<rttb::masks::boost::BoostMaskAccessor>(rtStructureSet->getStructure(j), geometricPtr);
 					CHECK_NO_THROW(boostMaskAccessorPtr->updateMask());
+
+					::boost::shared_ptr<core::GenericMaskedDoseIterator> spMaskedDoseIteratorTmp2 =
+					    ::boost::make_shared<core::GenericMaskedDoseIterator>(boostMaskAccessorPtr, doseAccessor1);
+					DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2);
+					rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(),
+					                               doseAccessor1->getDoseUID());
+					rttb::core::DVH dvh2 = *(calc2.generateDVH());
+
+					clock_t finish2(clock());
+					std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl;
+
+					CHECK_CLOSE(dvh.getMaximum(), dvh2.getMaximum(), 0.1);
+					CHECK_CLOSE(dvh.getMinimum(), dvh2.getMinimum(), 0.01);
+
 					//write the mask image to a file		
 					rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr);
 					CHECK(itkConverter2.process());
 					std::stringstream fileNameSstr2;
 					fileNameSstr2<<BoostMask_DIRNAME<<j<<".mhd";
 					std::cout << fileNameSstr2.str()<<std::endl;
 					rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage());
 					CHECK(writer2.writeITKFile());
 
-					clock_t finish2(clock());
-					std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl;
+					
 
 				}
 			}
 
 
 
 			RETURN_AND_REPORT_TEST_SUCCESS;
 
 
 		}
 
 	}//testing
 }//rttb