diff --git a/ShoulderCase/@DicomVolume/DicomVolume.m b/ShoulderCase/@DicomVolume/DicomVolume.m
index b4d1bf8..a756ade 100644
--- a/ShoulderCase/@DicomVolume/DicomVolume.m
+++ b/ShoulderCase/@DicomVolume/DicomVolume.m
@@ -1,55 +1,53 @@
 classdef (Abstract) DicomVolume < handle
 
   properties
     volume = []
     spatial = []
     dicomInfo = []
     dicomFolderPath = []
   end
 
 
 
   methods
     function set.volume(obj,volume)
       assert(length(size(volume))>=3,'Volume must be a 3D array')
-      obj.volume = squeeze(volume);
+      obj.volume = double(squeeze(volume));
     end
 
     function set.spatial(obj,spatial)
       requiredStructureField = {'PatientPositions','PixelSpacings','PatientOrientations'};
       assert(all(contains(fields(spatial),requiredStructureField)),'Invalid spatial structure');
       obj.spatial = spatial;
     end
 
     function loadDataFromFolder(obj,dicomFolderPath)
-      assert(obj.pathContainsDicomFiles(dicomFolderPath));
+      obj.assertPathContainsDicomFiles(dicomFolderPath);
       obj.dicomFolderPath = dicomFolderPath;
       [obj.volume,obj.spatial,~] = dicomreadVolume(dicomFolderPath);
       obj.loadDicomInfoFromFolder(dicomFolderPath);
     end
 
     function loadDicomInfoFromFolder(obj,dicomFolderPath)
-      assert(obj.pathContainsDicomFiles(dicomFolderPath));
+      obj.assertPathContainsDicomFiles(dicomFolderPath);
       dicomFiles = dir(fullfile(dicomFolderPath,'*.dcm'));
       obj.dicomInfo = dicominfo(fullfile(dicomFiles(1).folder,dicomFiles(1).name));
     end
 
-    function output = pathContainsDicomFiles(obj,path)
-      output = false;
+    function assertPathContainsDicomFiles(obj,path)
       assert(isfolder(path),'Provided argument is not a valid folder path');
       dicomFiles = dir(fullfile(path,'*.dcm'));
       assert(not(isempty(dicomFiles)),'No dicom file found there %s',path);
-      output = true;
     end
 
     function output = getPointIndexInVolume(obj,point)
       % the output are the three indices of the given point in the given volume.
       % i,j,k: output indices
       k = find(abs(obj.spatial.PatientPositions(:,3)-point(3)) == min(abs(obj.spatial.PatientPositions(:,3)-point(3))));
       i = round( (point(1) - obj.spatial.PatientPositions(k,1)) / obj.spatial.PixelSpacings(k,1) );
       j = round( (point(2) - obj.spatial.PatientPositions(k,2)) / obj.spatial.PixelSpacings(k,2) );
       assert(all([i j k] < size(obj.volume)),'Indices found are outside the volume boundaries');
       output = [i,j,k];
     end
   end
 end
diff --git a/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
index fd7f309..040aa57 100644
--- a/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
+++ b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
@@ -1,33 +1,36 @@
 classdef DicomVolumeNormaliser < handle & DicomVolume
 
   properties
     resolution = [];
     normalisedVolume = [];
   end
 
 
 
   methods
 
     function obj = DicomVolumeNormaliser(varargin)
       if (nargin == 1)
         obj.loadDataFromFolder(varargin{1});
       end
     end
 
     function setMinimalResolution(obj)
       zSpacing = obj.spatial.PatientPositions(2:end,3)-obj.spatial.PatientPositions(1:end-1,3);
       minZ = min(zSpacing,[],'all');
       minXY = min(obj.spatial.PixelSpacings,[],'all');
       obj.resolution = min([minXY minZ],[],'all')/2;
     end
 
     function normalise(obj)
       normalisedSizeZ = round(abs(obj.spatial.PatientPositions(1,3) - obj.spatial.PatientPositions(end,3))/obj.resolution);
       normalisedSizeX = round(max(obj.spatial.PixelSpacings(:,1),[],'all')*size(obj.volume,1)/obj.resolution);
       normalisedSizeY = round(max(obj.spatial.PixelSpacings(:,2),[],'all')*size(obj.volume,2)/obj.resolution);
       obj.normalisedVolume = imresize3(obj.volume,[normalisedSizeX normalisedSizeY normalisedSizeZ]);
+      newZPositions = (obj.spatial.PatientPositions(1,3):obj.resolution:obj.spatial.PatientPositions(end,3))';
+      obj.spatial.PatientPositions = [repmat(obj.spatial.PatientPositions(1,1:2),length(newZPositions),1) newZPositions];
+      obj.spatial.PixelSpacings = repmat(obj.resolution,size(obj.spatial.PatientPositions,1),2);
     end
   end
 
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
index bee14b4..2e9cfa0 100644
--- a/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
+++ b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
@@ -1,30 +1,90 @@
 classdef DicomVolumeSlicer < handle & DicomVolume & DicomVolumeNormaliser
 
   properties
     sliced
-    slicedPixelSpacing
   end
 
-
+  properties (Access = public)
+    slicedPixelSpacings
+    slicedPlane
+    slicedX
+    slicedY
+    slicedZ
+  end
 
   methods
 
     function obj = DicomVolumeSlicer(varargin)
       if (nargin == 1)
         obj.loadDataFromFolder(varargin{1});
       end
     end
 
     function slice(obj,point,normalVector)
+      obj.slicedPlane = Plane();
+      obj.slicedPlane.point = point;
+      obj.slicedPlane.normal = normalVector;
       pointIndices = obj.getPointIndexInVolume(point);
-      vectorIndices = pointIndices - obj.getPointIndexInVolume(point+normalVector);
-
+      vectorIndices = obj.getPointIndexInVolume(point+50*normalVector/norm(normalVector)) - pointIndices;
       [obj.sliced,x,y,z] = obliqueslice(obj.volume,pointIndices,vectorIndices,...
-                                'FillValues',min(obj.volume,[],'all'));
-      obj.sliced = obj.sliced';
-      obj.sliced = obj.sliced*obj.dicomInfo.RescaleSlope + obj.dicomInfo.RescaleIntercept;
+                                        'FillValues',-2000);
+      obj.slicedX = x;
+      obj.slicedY = y;
+      obj.slicedZ = z;
+      obj.measureSlicedPixelSpacings;
+      obj.sliced(obj.sliced <= -1999) = nan;
+    end
+
+    function [vectorIndices] = orientSliceUpwardVector(obj,point,vector)
+      pointIndices = obj.getPointIndexInSlice(point);
+      vectorIndices = obj.getPointIndexInSlice(obj.slicedPlane.projectOnPlane(point + 10*vector/norm(vector))) - pointIndices;
+      rotation = vrrotvec([vectorIndices 0],[-1 0 0]);
+      obj.sliced = imrotate(obj.sliced,rotation(3)*180*rotation(4)/pi);
+      obj.slicedX = imrotate(obj.slicedX,rotation(3)*180*rotation(4)/pi);
+      obj.slicedY = imrotate(obj.slicedY,rotation(3)*180*rotation(4)/pi);
+      obj.slicedZ = imrotate(obj.slicedZ,rotation(3)*180*rotation(4)/pi);
+      obj.measureSlicedPixelSpacings;
+    end
+
+    function [top,bottom,left,right] = crop(obj,center,width,height)
+      % center: point in volume
+      % width, height: length in cm
+      centerIndex = obj.getPointIndexInSlice(center);
+
+      top = centerIndex(1)-round((height/2)/obj.slicedPixelSpacings(2));
+      bottom = centerIndex(1)+round((height/2)/obj.slicedPixelSpacings(2));
+      left = centerIndex(1)-round((width/2)/obj.slicedPixelSpacings(1));
+      right = centerIndex(1)+round((width/2)/obj.slicedPixelSpacings(1));
+
+      % obj.sliced = obj.sliced(top:bottom,left:right);
+      % obj.slicedX = obj.slicedX(top:bottom,left:right);
+      % obj.slicedY = obj.slicedY(top:bottom,left:right);
+      % obj.slicedZ = obj.slicedZ(top:bottom,left:right);
+    end
+
+    function output = getPointIndexInSlice(obj,point)
+      pointIndices = obj.getPointIndexInVolume(point);
+      [i,j] = ind2sub(size(obj.sliced),find(round(obj.slicedX) == pointIndices(1) &...
+                                            round(obj.slicedY) == pointIndices(2) &...
+                                            round(obj.slicedZ) == pointIndices(3)));
+      assert(not(isempty([i,j])),'Point has not been found in slice');
+      output = [i(1),j(1)];
+    end
+
+    function [originIndex, bottomLeftIndex, topRightIndex] = measureSlicedPixelSpacings(obj)
+      originIndex = [round(obj.slicedX(1,1)),...
+                     round(obj.slicedY(1,1)),...
+                     round(obj.slicedZ(1,1))];
+      bottomLeftIndex = [round(obj.slicedX(end,1)),...
+                         round(obj.slicedY(end,1)),...
+                         round(obj.slicedZ(end,1))];
+      topRightIndex = [round(obj.slicedX(1,end)),...
+                       round(obj.slicedY(1,end)),...
+                       round(obj.slicedZ(1,end))];
+      % obj.slicedPixelSpacings = [norm(obj.volume(bottomLeftIndex)-obj.volume(originIndex))/size(obj.sliced,1),...
+      %                            norm(obj.volume(topRightIndex)-obj.volume(originIndex))/size(obj.sliced,2)];
     end
 
   end
 
 end