diff --git a/utils/BinaryVolume/@BinaryCylinder/BinaryCylinder.m b/utils/BinaryVolume/@BinaryCylinder/BinaryCylinder.m
new file mode 100644
index 0000000..9bd34cc
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryCylinder/BinaryCylinder.m
@@ -0,0 +1,28 @@
+classdef BinaryCylinder < BinaryVolume
+
+    methods (Access = public)
+
+        function obj = BinaryCylinder(radius,length,varargin)
+            if nargin == 3
+                center = varargin{1};
+            else
+                center = [0 0 0];
+            end
+            obj@BinaryVolume([2*radius 2*radius length]);
+            obj.translate(center);
+
+            squaredDistanceToZAxis = obj.getSquaredDistanceToPointXY(center);
+            cylinder = squaredDistanceToZAxis <= radius^2;
+
+            cylinderStartPoint = obj.getIndicesOfPointInVolume(center - [0 0 length/2]);
+            cylinderEndPoint = obj.getIndicesOfPointInVolume(center + [0 0 length/2]);
+            cylinder(:,:,1:cylinderStartPoint(3)) = false;
+            cylinder(:,:,cylinderEndPoint(3):end) = false;
+            obj.volume = cylinder;
+
+            obj.resize(obj.size + 2*obj.resolution);
+        end
+
+    end
+
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryDicomVolume/BinaryDicomVolume.m b/utils/BinaryVolume/@BinaryDicomVolume/BinaryDicomVolume.m
new file mode 100644
index 0000000..54226d6
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryDicomVolume/BinaryDicomVolume.m
@@ -0,0 +1,38 @@
+classdef BinaryDicomVolume < BinaryVolume
+
+  properties (Access = public)
+    dicomVolume
+  end
+
+
+  methods 
+
+    function obj = BinaryDicomVolume(dicomFolderPath)
+      [volume,spatial,~] = dicomreadVolume(dicomFolderPath);
+      volume = squeeze(volume);
+      resolution = getSmallestResolutionDividedBy2(spatial);
+      volume = normaliseDicomVolume(volume,resolution);
+      
+
+      obj@BinaryVolume(size(volume));
+      obj.dicomVolume = volume;
+      obj.setResolution(resolution);
+    end
+
+  end
+
+end
+
+
+function output = getSmallestResolutionDividedBy2(spatial)
+  resolution(1) = min(spatial.PixelSpacings(:,1),[],'all')/2;
+  resolution(2) = min(spatial.PixelSpacings(:,2),[],'all')/2;
+  resolution(3) = (spatial.PatientPositions(end,3) - spatial.PatientPositions(1,3))/(2*size(spatial.PatientPositions,1));
+  output = resolution;
+end
+
+function output = normaliseDicomVolume(volume,newResolution)
+  newSize = abs(size(volume)./newResolution);
+  volume = imresize3(volume,newSize);
+  output = volume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryDicomVolume/filterHUBelow.m b/utils/BinaryVolume/@BinaryDicomVolume/filterHUBelow.m
new file mode 100644
index 0000000..d55b630
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryDicomVolume/filterHUBelow.m
@@ -0,0 +1,4 @@
+function output = filterHUBelow(obj,threshold)
+  obj.setVolume(obj.dicomVolume < threshold);
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryDicomVolume/filterHUBetween.m b/utils/BinaryVolume/@BinaryDicomVolume/filterHUBetween.m
new file mode 100644
index 0000000..a30433a
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryDicomVolume/filterHUBetween.m
@@ -0,0 +1,3 @@
+function output = filterHUBetween(obj,threshold)
+  obj.setVolume( (obj.dicomVolume > threshold(1)) & (obj.dicomVolume < threshold(2)) );
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryDicomVolume/filterHUOver.m b/utils/BinaryVolume/@BinaryDicomVolume/filterHUOver.m
new file mode 100644
index 0000000..bbc85b4
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryDicomVolume/filterHUOver.m
@@ -0,0 +1,4 @@
+function output = filterHUOver(obj,threshold)
+  obj.setVolume(obj.dicomVolume > threshold);
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryRectangle/BinaryRectangle.m b/utils/BinaryVolume/@BinaryRectangle/BinaryRectangle.m
new file mode 100644
index 0000000..1a30372
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryRectangle/BinaryRectangle.m
@@ -0,0 +1,20 @@
+classdef BinaryRectangle < BinaryVolume
+
+    methods (Access = public)
+
+        function obj = BinaryRectangle(rectangleSize,varargin)
+            if nargin == 2
+                center = varargin{1};
+            else
+                center =  [0 0 0];
+            end
+            obj@BinaryVolume(rectangleSize)
+            obj.translate(center);
+            obj.fillVolume();
+            
+            obj.resize(obj.size + 2*obj.resolution);
+        end
+        
+    end
+
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinarySphere/BinarySphere.m b/utils/BinaryVolume/@BinarySphere/BinarySphere.m
new file mode 100644
index 0000000..4074b17
--- /dev/null
+++ b/utils/BinaryVolume/@BinarySphere/BinarySphere.m
@@ -0,0 +1,21 @@
+classdef BinarySphere < BinaryVolume
+
+    methods (Access = public)
+
+        function obj = BinarySphere(radius,varargin)
+            if nargin == 2
+                center = varargin{1};
+            else
+                center = [0 0 0];
+            end
+            obj@BinaryVolume(repmat(2*radius,[1 3]));
+            obj.translate(center);
+            squaredDistanceToCenter = obj.getSquaredDistanceToPointXYZ(center);
+            obj.volume = squaredDistanceToCenter <= radius^2;
+            
+            obj.resize(obj.size + 2*obj.resolution);
+        end
+
+    end
+
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/BinaryVolume.m b/utils/BinaryVolume/@BinaryVolume/BinaryVolume.m
new file mode 100644
index 0000000..df1639b
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/BinaryVolume.m
@@ -0,0 +1,35 @@
+classdef BinaryVolume < handle
+
+  properties (Access = public)
+    % origin is the spatial position of volume(1,1,1)
+    origin = [-50 -50 -50];
+    size = [100 100 100];
+    resolution = [1 1 1];
+    volume = false(100,100,100);
+  end
+
+
+
+  methods 
+
+    function obj = BinaryVolume(varargin)
+      if nargin == 1
+        obj.volume = false(varargin{1});
+        obj.origin = -obj.size/2;
+      end
+    end
+
+    function set.volume(obj,value)
+      obj.volume = value;
+      obj.size = size(obj.volume).*obj.resolution;
+    end
+
+    function set.resolution(obj,value)
+      assert(all(value > 0));
+      obj.resolution = value;
+      obj.size = size(obj.volume).*obj.resolution;
+    end
+
+  end
+
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/addBoundingBox.m b/utils/BinaryVolume/@BinaryVolume/addBoundingBox.m
new file mode 100644
index 0000000..1c92ce4
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/addBoundingBox.m
@@ -0,0 +1,18 @@
+function output = addBoundingBox(obj)
+  obj.volume(:,1,1) = true;
+  obj.volume(:,end,1) = true;
+  obj.volume(:,1,end) = true;
+  obj.volume(:,end,end) = true;
+
+  obj.volume(1,:,1) = true;
+  obj.volume(end,:,1) = true;
+  obj.volume(1,:,end) = true;
+  obj.volume(end,:,end) = true;
+
+  obj.volume(1,1,:) = true;
+  obj.volume(end,1,:) = true;
+  obj.volume(1,end,:) = true;
+  obj.volume(end,end,:) = true;
+
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/changeResolution.m b/utils/BinaryVolume/@BinaryVolume/changeResolution.m
new file mode 100644
index 0000000..3a0245b
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/changeResolution.m
@@ -0,0 +1,10 @@
+function output = changeResolution(obj,newResolution)
+  assert(all(newResolution > 0));
+  if all(newResolution == obj.resolution)
+    output = obj;
+    return
+  end
+  obj.volume = imresize3(obj.volume,ceil(obj.size./newResolution));
+  obj.resolution = newResolution;
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/clearVolume.m b/utils/BinaryVolume/@BinaryVolume/clearVolume.m
new file mode 100644
index 0000000..6d41d75
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/clearVolume.m
@@ -0,0 +1,4 @@
+function output = clearVolume(obj)
+    obj.volume = false(size(obj.volume));
+    output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/containsVolume.m b/utils/BinaryVolume/@BinaryVolume/containsVolume.m
new file mode 100644
index 0000000..028e45c
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/containsVolume.m
@@ -0,0 +1,17 @@
+function output = containsVolume(obj,binaryVolume)
+  boundaryMin = obj.origin;
+  boundaryMax = obj.origin + obj.size;
+  testedVolumeMin = binaryVolume.origin;
+  testedVolumeMax = binaryVolume.origin + binaryVolume.size;
+
+  if not( all(boundaryMin <= testedVolumeMin) )
+    output = false;
+    return
+  end
+  if not( all(boundaryMax >= testedVolumeMax) )
+    output = false;
+    return
+  end
+
+  output = true;
+end
diff --git a/utils/BinaryVolume/@BinaryVolume/copy.m b/utils/BinaryVolume/@BinaryVolume/copy.m
new file mode 100644
index 0000000..e5beef6
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/copy.m
@@ -0,0 +1,8 @@
+function output = copy(obj)
+  copiedVolume = BinaryVolume();
+  copiedVolume.setResolution(obj.resolution);
+  copiedVolume.setVolume(obj.volume);
+  copiedVolume.translate(obj.origin-copiedVolume.origin);
+
+  output = copiedVolume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/difference.m b/utils/BinaryVolume/@BinaryVolume/difference.m
new file mode 100644
index 0000000..5918a0d
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/difference.m
@@ -0,0 +1,4 @@
+function output = difference(obj,volumeToSubstract)
+    output = copy(obj);
+    output.volume = output.volume & not(volumeToSubstract.volume);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/exportPly.m b/utils/BinaryVolume/@BinaryVolume/exportPly.m
new file mode 100644
index 0000000..99e76b6
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/exportPly.m
@@ -0,0 +1,9 @@
+function exportPly(obj,filename)
+  volumeToExport = bwmorph3(obj.volume,'remove');
+  [X,Y,Z] = ind2sub(size(volumeToExport),find(volumeToExport));
+  X = X*obj.resolution(1) + obj.origin(1);
+  Y = Y*obj.resolution(2) + obj.origin(2);
+  Z = Z*obj.resolution(3) + obj.origin(3);
+  points = ([X Y Z]);
+  plywrite(filename,points);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/fillVolume.m b/utils/BinaryVolume/@BinaryVolume/fillVolume.m
new file mode 100644
index 0000000..283a4fb
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/fillVolume.m
@@ -0,0 +1,4 @@
+function output = fillVolume(obj)
+    obj.volume = true(size(obj.volume));
+    output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/getIndicesOfPointInVolume.m b/utils/BinaryVolume/@BinaryVolume/getIndicesOfPointInVolume.m
new file mode 100644
index 0000000..f5bccf6
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/getIndicesOfPointInVolume.m
@@ -0,0 +1,3 @@
+function output = getIndicesOfPointInVolume(obj,point)
+  output = round((point - obj.origin)./obj.resolution) + [1 1 1];
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/getSmallestCommonBoundingVolume.m b/utils/BinaryVolume/@BinaryVolume/getSmallestCommonBoundingVolume.m
new file mode 100644
index 0000000..75659c4
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/getSmallestCommonBoundingVolume.m
@@ -0,0 +1,7 @@
+function output = getSmallestCommonBoundingVolume(obj,comparedVolume)
+  commonOrigin = min(obj.origin,comparedVolume.origin);
+  commonDiagonal = max( (obj.origin + obj.size) , (comparedVolume.origin + comparedVolume.size) );
+  boundingVolume = BinaryVolume(ceil(commonDiagonal-commonOrigin));
+  boundingVolume.setOrigin(commonOrigin);
+  output = boundingVolume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXY.m b/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXY.m
new file mode 100644
index 0000000..77d03fd
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXY.m
@@ -0,0 +1,10 @@
+function output = getSquaredDistanceToPointXY(obj,point)
+  [X,Y] = ndgrid(int32(1:size(obj.volume,1)),int32(1:size(obj.volume,2)),int32(1:size(obj.volume,3)));
+
+  point = obj.getIndicesOfPointInVolume(point);
+
+  X = (X - point(1)) * obj.resolution(1);
+  Y = (Y - point(2)) * obj.resolution(2);
+
+  output = X.^2 + Y.^2;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXYZ.m b/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXYZ.m
new file mode 100644
index 0000000..28915b1
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/getSquaredDistanceToPointXYZ.m
@@ -0,0 +1,11 @@
+function output = getSquaredDistanceToPointXYZ(obj,point)
+  [X,Y,Z] = ndgrid(int32(1:size(obj.volume,1)),int32(1:size(obj.volume,2)),int32(1:size(obj.volume,3)));
+
+  point = obj.getIndicesOfPointInVolume(point);
+
+  X = (X - point(1)) * obj.resolution(1);
+  Y = (Y - point(2)) * obj.resolution(2);
+  Z = (Z - point(3)) * obj.resolution(3);
+
+  output = X.^2 + Y.^2 + Z.^2;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/getTriangulation.m b/utils/BinaryVolume/@BinaryVolume/getTriangulation.m
new file mode 100644
index 0000000..f034067
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/getTriangulation.m
@@ -0,0 +1,5 @@
+function output = getTriangulation(obj)
+    volume = bwmorph3(obj.volume,'remove');
+    [X,Y,Z] = ind2sub(size(obj.volume),find(volume));
+    output = delaunayTriangulation(X,Y,Z);
+end
diff --git a/utils/BinaryVolume/@BinaryVolume/insertVolume.m b/utils/BinaryVolume/@BinaryVolume/insertVolume.m
new file mode 100644
index 0000000..ea36cd4
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/insertVolume.m
@@ -0,0 +1,32 @@
+function output = insertVolume(obj,insertedVolume)
+  copiedInsertedVolume = insertedVolume.copy();
+  startingResolution = obj.resolution;
+  lowestResolution = min(obj.resolution,copiedInsertedVolume.resolution);
+  obj.changeResolution(lowestResolution);
+  copiedInsertedVolume.changeResolution(lowestResolution);
+
+  insertionPoint = obj.getIndicesOfPointInVolume(copiedInsertedVolume.origin);
+  insertionPointOpposite = (insertionPoint + size(copiedInsertedVolume.volume) - [1 1 1]);
+  insertionPoint = min(size(obj.volume),max([1 1 1],insertionPoint));
+  insertionPointOpposite = min(size(obj.volume),max([1 1 1],insertionPointOpposite));
+  if any(insertionPoint == insertionPointOpposite)
+    return
+  end
+  extractionPoint = copiedInsertedVolume.getIndicesOfPointInVolume(obj.origin + obj.resolution.*(insertionPoint-[1 1 1]));
+  extractionPointOpposite = extractionPoint + (insertionPointOpposite - insertionPoint);
+  extractionPoint = min(size(copiedInsertedVolume.volume),max([1 1 1],extractionPoint));
+  extractionPointOpposite = min(size(copiedInsertedVolume.volume),max([1 1 1],extractionPointOpposite));
+  
+  insertionPointOpposite = insertionPoint + (extractionPointOpposite - extractionPoint);
+
+
+  subvolumeToInsert = copiedInsertedVolume.volume(extractionPoint(1):extractionPointOpposite(1),...
+    extractionPoint(2):extractionPointOpposite(2),...
+    extractionPoint(3):extractionPointOpposite(3)); 
+  obj.volume(insertionPoint(1):insertionPointOpposite(1),...
+    insertionPoint(2):insertionPointOpposite(2),...
+    insertionPoint(3):insertionPointOpposite(3)) = subvolumeToInsert;
+
+  obj.changeResolution(startingResolution);
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/intersection.m b/utils/BinaryVolume/@BinaryVolume/intersection.m
new file mode 100644
index 0000000..9e95b56
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/intersection.m
@@ -0,0 +1,4 @@
+function output = intersection(obj,volumeToIntersect)
+    output = copy(obj);
+    output.volume = output.volume & volumeToIntersect.volume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/isempty.m b/utils/BinaryVolume/@BinaryVolume/isempty.m
new file mode 100644
index 0000000..c2bcdff
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/isempty.m
@@ -0,0 +1,3 @@
+function output = isempty(obj)
+    output = numel(find(obj.volume)) == 0;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/minus.m b/utils/BinaryVolume/@BinaryVolume/minus.m
new file mode 100644
index 0000000..3edc4f5
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/minus.m
@@ -0,0 +1,6 @@
+function output = minus(obj,volumeSubstracted)
+  currentVolumeForDifference = insertVolume(getSmallestCommonBoundingVolume(obj,volumeSubstracted),obj);
+  substractedVolumeForDifference = insertVolume(getSmallestCommonBoundingVolume(obj,volumeSubstracted),volumeSubstracted);
+
+  output = difference(currentVolumeForDifference,substractedVolumeForDifference);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/mtimes.m b/utils/BinaryVolume/@BinaryVolume/mtimes.m
new file mode 100644
index 0000000..a08393b
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/mtimes.m
@@ -0,0 +1,6 @@
+function output = mtimes(obj,volumeCompared)
+  currentVolumeForIntersection = insertVolume(getSmallestCommonBoundingVolume(obj,volumeCompared),obj);
+  comparedVolumeForIntersection = insertVolume(getSmallestCommonBoundingVolume(obj,volumeCompared),volumeCompared);
+
+  output = intersection(currentVolumeForIntersection,comparedVolumeForIntersection);
+end
diff --git a/utils/BinaryVolume/@BinaryVolume/plus.m b/utils/BinaryVolume/@BinaryVolume/plus.m
new file mode 100644
index 0000000..50ec09a
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/plus.m
@@ -0,0 +1,6 @@
+function output = plus(obj,volumeAdded)
+  currentVolumeForUnion = insertVolume(getSmallestCommonBoundingVolume(obj,volumeAdded),obj);
+  addedVolumeForUnion = insertVolume(getSmallestCommonBoundingVolume(obj,volumeAdded),volumeAdded);
+
+  output = union(currentVolumeForUnion,addedVolumeForUnion);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/removeEmptySpace.m b/utils/BinaryVolume/@BinaryVolume/removeEmptySpace.m
new file mode 100644
index 0000000..6dcbe13
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/removeEmptySpace.m
@@ -0,0 +1,19 @@
+function output = removeEmptySpace(obj)
+    [X,Y,Z] = ind2sub(size(obj.volume),find(obj.volume));
+    if isempty([X,Y,Z])
+        lowestPoint = [1 1 1];
+        opposite = [2 2 2];
+    else
+        lowestPoint = min([X Y Z])
+        opposite = max([X,Y,Z])
+    end
+
+    resizedVolume = BinaryVolume([opposite-lowestPoint]);
+    resizedVolume.setOrigin(obj.origin + lowestPoint.*obj.resolution);
+    resizedVolume.insertVolume(obj);
+
+    obj.volume = resizedVolume.volume;
+    obj.origin = resizedVolume.origin;
+
+    output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/resize.m b/utils/BinaryVolume/@BinaryVolume/resize.m
new file mode 100644
index 0000000..c887ff6
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/resize.m
@@ -0,0 +1,12 @@
+function output = resize(obj,newSize)
+  currentCenter = obj.origin + obj.size/2;
+  resizedEmptyVolume = BinaryVolume(newSize);
+  resizedEmptyVolume.translate(currentCenter);
+  resizedVolume = resizedEmptyVolume.insertVolume(obj);
+
+  obj.setVolume(resizedVolume.volume);
+  obj.setOrigin(resizedVolume.origin);
+  obj.setResolution(resizedVolume.resolution);
+
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/rotate.m b/utils/BinaryVolume/@BinaryVolume/rotate.m
new file mode 100644
index 0000000..b476059
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/rotate.m
@@ -0,0 +1,8 @@
+function output = rotate(obj,rotationVector)
+  rotationVector(4) = rotationVector(4)*180/pi;
+  startingResolution = obj.resolution;
+  obj.setResolution([1 1 1]);
+  obj.volume = imrotate3(obj.volume,rotationVector(4),rotationVector(1:3),'nearest','loose');
+  obj.setResolution(startingResolution);
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/setOrigin.m b/utils/BinaryVolume/@BinaryVolume/setOrigin.m
new file mode 100644
index 0000000..d955d2d
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/setOrigin.m
@@ -0,0 +1,3 @@
+function setOrigin(obj,newOrigin)
+  obj.origin = newOrigin;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/setResolution.m b/utils/BinaryVolume/@BinaryVolume/setResolution.m
new file mode 100644
index 0000000..8e67fbd
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/setResolution.m
@@ -0,0 +1,3 @@
+function setResolution(obj,newResolution)
+  obj.resolution = newResolution;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/setVolume.m b/utils/BinaryVolume/@BinaryVolume/setVolume.m
new file mode 100644
index 0000000..a541898
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/setVolume.m
@@ -0,0 +1,3 @@
+function setVolume(obj,volume)
+  obj.volume = volume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/show.m b/utils/BinaryVolume/@BinaryVolume/show.m
new file mode 100644
index 0000000..edaf622
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/show.m
@@ -0,0 +1,3 @@
+function show(obj)
+  volshow(obj.volume);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/translate.m b/utils/BinaryVolume/@BinaryVolume/translate.m
new file mode 100644
index 0000000..510711c
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/translate.m
@@ -0,0 +1,4 @@
+function output = translate(obj,translation)
+  obj.origin = obj.origin + translation;
+  output = obj;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/uminus.m b/utils/BinaryVolume/@BinaryVolume/uminus.m
new file mode 100644
index 0000000..6538556
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/uminus.m
@@ -0,0 +1,5 @@
+function output = uminus(obj)
+    output = copy(obj);
+    output.volume = not(output.volume);
+    output.resize(output.size + 2*output.resolution);
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/@BinaryVolume/union.m b/utils/BinaryVolume/@BinaryVolume/union.m
new file mode 100644
index 0000000..80173ab
--- /dev/null
+++ b/utils/BinaryVolume/@BinaryVolume/union.m
@@ -0,0 +1,4 @@
+function output = union(obj,volumeToUnite)
+    output = copy(obj);
+    output.volume = output.volume | volumeToUnite.volume;
+end
\ No newline at end of file
diff --git a/utils/BinaryVolume/plywrite.m b/utils/BinaryVolume/plywrite.m
new file mode 100644
index 0000000..860a576
--- /dev/null
+++ b/utils/BinaryVolume/plywrite.m
@@ -0,0 +1,74 @@
+function plywrite(filename,points,varargin)
+  fid = fopen(filename,'w');
+
+  % pointsTriangulation = delaunayTriangulation(points);
+  % points = pointsTriangulation.Points;
+  % facets = getFacetsFromTetrahedrons(pointsTriangulation.ConnectivityList);
+  % facets = getFacetsFromTetrahedrons(alphaTriangulation(alphaShape(points),1));
+  facets = boundary(points,1) - [1 1 1];
+
+  writeHeader(fid,size(points,1),size(facets,1));
+  if nargin == 2
+    writePoints(fid,points);
+  else
+    colormap = varargin{1};
+    writePointsWithColormap(fid,points,colormap);
+  end
+  writeFacets(fid,facets);
+  fclose(fid);
+  disp('Ply successfully written.');
+end
+
+
+function output = getFacetsFromTetrahedrons(connectivity)
+  connectivity = connectivity - [1 1 1 1];
+  facets = [connectivity(:,[1 2 3]);connectivity(:,[1 2 4]);connectivity(:,[1 3 4]);connectivity(:,[2 3 4])];
+  facets = sort(facets,2);
+  facets = sortrows(facets);
+  duplicates = find(all(diff(facets) == 0,2));
+  facets([duplicates;duplicates + 1],:) = [];
+  output = facets;
+end
+
+
+function writeHeader(fid,numberOfPoints,numberOfFacets)
+  fprintf(fid,...
+    ['ply\n'...
+    'format ascii 1.0\n'...
+    'element vertex %d\n'...
+    'property float32 x\n'...
+    'property float32 y\n'...
+    'property float32 z\n'...
+    'element face %d\n'...
+    'property list uint8 int32 vertex_indices\n'...
+    'end_header\n'],...
+    numberOfPoints,numberOfFacets);
+end
+
+
+function writePoints(fid,points)
+  for i = 1:size(points,1)
+    fprintf(fid,...
+      '%.6f %.6f %.6f\n',...
+      points(i,1),points(i,2),points(i,3));
+  end
+end
+
+
+function writePointsWithColormap(fid,points,colormap)
+  for i = 1:size(points,1)
+    fprintf(fid,...
+      '%.6f %.6f %.6f %u %u %u\n',...
+      points(i,1),points(i,2),points(i,3),...
+      colormap(i,1),colormap(i,2),colormap(i,3));
+  end
+end
+
+
+function writeFacets(fid,facets)
+  for i = 1:size(facets,1)
+    fprintf(fid,...
+      '%d %d %d %d\n',...
+      size(facets,2),facets(i,1),facets(i,2),facets(i,3));
+  end
+end
\ No newline at end of file