diff --git a/ShoulderCase/@Muscles/getDicomPath.m b/ShoulderCase/@Muscles/getDicomPath.m
new file mode 100644
index 0000000..d62c379
--- /dev/null
+++ b/ShoulderCase/@Muscles/getDicomPath.m
@@ -0,0 +1,16 @@
+function output = getDicomPath(obj)
+  output = [];
+  SCase = obj.shoulder.SCase;
+  rawDicomPath = fullfile(SCase.dataCTPath,'dicom');
+  smoothDicomPath = fullfile([SCase.dataCTPath(1:end-1) '2'],'dicom');
+
+  if exist(smoothDicomPath,'dir')
+    output = smoothDicomPath;
+  elseif exist(rawDicomPath,'dir')
+    output = rawDicomPath;
+  else
+    raisedError.identifier = 'ShoulderCase:Shoulder:Muscles:getDicomPath:InvalidPath';
+    raisedError.message = '%s is not a valid folder';
+    error(raisedError);
+  end
+end
diff --git a/ShoulderCase/@Muscles/getRotatorCuffSlice.m b/ShoulderCase/@Muscles/getRotatorCuffSlice.m
new file mode 100644
index 0000000..63f4ad1
--- /dev/null
+++ b/ShoulderCase/@Muscles/getRotatorCuffSlice.m
@@ -0,0 +1,13 @@
+function output = getRotatorCuffSlice(obj)
+  output = [];
+  dicomPath = obj.getDicomPath;
+  dicomFiles = dir(dicomPath);
+  dicomInfo = dicominfo(dicomFiles(3).name);
+  [volume,spatial,dimension] = dicomreadVolume(dicomPath);
+
+  volume = squeeze(volume);
+  scapulaOriginIndices = getPointIndexInDicomVolume(obj.shoulder.scapula.coordSys.origin,volume,spatial);
+  slice = obliqueslice(volume,scapulaOriginIndices,obj.shoulder.scapula.coordSys.ML)
+
+  output = slice;
+end
diff --git a/ShoulderCase/getPointIndexInDicomVolume.m b/ShoulderCase/getPointIndexInDicomVolume.m
new file mode 100644
index 0000000..a4da26c
--- /dev/null
+++ b/ShoulderCase/getPointIndexInDicomVolume.m
@@ -0,0 +1,8 @@
+function output = getPointIndexInDicomVolume(point,volume,spatial)
+  % the output are the three indices of the given point in the given volume.
+  % i,j,k: output indices
+  k = find(abs(spatial.PatientPositions(:,3)-point(3)) == min(abs(spatial.PatientPositions(:,3)-point(3))));
+  i = round( (point(1) - spatial.PatientPositions(k,1)) / spatial.PixelSpacings(k,1) );
+  j = round( (point(2) - spatial.PatientPositions(k,2)) / spatial.PixelSpacings(k,2) );
+  output = [i,j,k];
+end