diff --git a/ShoulderCase/@MaskSubdivider/MaskSubdivider.m b/ShoulderCase/@MaskSubdivider/MaskSubdivider.m
new file mode 100644
index 0000000..2fded4f
--- /dev/null
+++ b/ShoulderCase/@MaskSubdivider/MaskSubdivider.m
@@ -0,0 +1,119 @@
+classdef MaskSubdivider < handle
+
+  properties (Access = private)
+    % 10° slices
+    numberOfAngularDivisions = 36;
+    numberOfRadialDivisions = 10;
+    center
+    mask
+
+    sinusGrid
+    cosineGrid
+    angleGrid
+    radiusGrid
+  end
+
+  methods 
+    function obj = MaskSubdivider(mask)
+      obj.mask = mask;
+      % default center is the middle of the mask
+      obj.center = round(size(obj.mask)/2);
+    end
+
+    function setNumberOfAngularDivisions(obj, value)
+      obj.numberOfAngularDivisions = value;
+    end
+
+    function setNumberOfRadialDivisions(obj, value)
+      obj.numberOfRadialDivisions = value;
+    end
+
+    function setCenter(obj, centerIndices)
+      obj.center = centerIndices;
+    end
+
+    function set.center(obj, value)
+      obj.center = value;
+      obj.updateGrids();
+    end
+
+    function updateGrids(obj)
+      [xGrid, yGrid] = meshgrid(1:size(obj.mask, 1), 1:size(obj.mask, 2));
+      xGrid = xGrid-obj.center(1);
+      yGrid = yGrid-obj.center(2);
+      obj.radiusGrid = sqrt(xGrid.^2 + yGrid.^2);
+      obj.sinusGrid = yGrid./obj.radiusGrid;
+      obj.cosineGrid = xGrid./obj.radiusGrid;
+      partialAngleGrid = acosd(obj.cosineGrid);
+      obj.angleGrid = [...
+        (360 - partialAngleGrid).*(obj.sinusGrid < 0)...
+        + partialAngleGrid.*(obj.sinusGrid >= 0)];
+    end
+
+    function output = getAngularDivision(obj, middleAngle, angularSpan)
+      minAngle = mod(middleAngle-angularSpan/2, 360);
+      maxAngle = mod(middleAngle+angularSpan/2, 360);
+      output = obj.angleGrid > minAngle & obj.angleGrid < maxAngle;
+    end
+
+    function output = getRadialDivision(obj, middleRadius, radiusSpan)
+      minRadius = middleRadius - radiusSpan/2;
+      maxRadius = middleRadius + radiusSpan/2;
+      output = obj.radiusGrid > minRadius & obj.radiusGrid < maxRadius;
+    end
+
+    function output = getMaskSubdivisions(obj)
+      maskSubdivisions = [];
+
+      [minRadius, maxRadius] = obj.getMinMaxGridValueInMask(obj.radiusGrid);
+      radialIncrement = (maxRadius - minRadius)/obj.numberOfRadialDivisions;
+      radialIteratorStart = minRadius + radialIncrement/2;
+      radialIteratorEnd = maxRadius - radialIncrement/2;
+      for radius = radialIteratorStart:radialIncrement:radialIteratorEnd
+        radialDivision = obj.getRadialDivision(radius, radialIncrement);
+
+        [minAngle, maxAngle] = obj.getMinMaxGridValueInMask(obj.angleGrid);
+        angularIncrement = (maxAngle - minAngle)/obj.numberOfAngularDivisions;
+        angularIteratorStart = minAngle + angularIncrement/2;
+        angularIteratorEnd = maxAngle - angularIncrement/2;
+        for angle = angularIteratorStart:angularIncrement:angularIteratorEnd
+          angularDivision = obj.getAngularDivision(angle, angularIncrement);
+          maskSubdivisions(:,:,end+1) =...
+            obj.mask & radialDivision & angularDivision;
+        end
+      end
+      output = maskSubdivisions;
+      output = obj.removeEmptySubdivisions(maskSubdivisions);
+    end
+
+    function [minOutput, maxOutput] = getMinMaxGridValueInMask(obj, inputGrid)
+      gridInMask = obj.mask.*inputGrid;
+      gridInMask(not(obj.mask)) = nan;
+      minOutput = min(gridInMask, [], "all");
+      maxOutput = max(gridInMask, [], "all");
+    end
+
+    function output = removeEmptySubdivisions(obj, subdivisionsArray)
+      subdivisionsToRemove = [];
+      for i = 1:size(subdivisionsArray, 3)
+        if not(any(subdivisionsArray(:,:,i), "all"))
+          subdivisionsToRemove(end+1) = i;
+        end
+      end
+      subdivisionsArray(:,:,subdivisionsToRemove) = [];
+      output = subdivisionsArray;
+    end
+
+    function output = getMaskSubdivisionsStacked(obj)
+      maskSubdivisions = obj.getMaskSubdivisions();
+      maskSubdivisionsStacked = false(size(obj.mask));
+      for i = 1:size(maskSubdivisions, 3)
+        maskSubdivisionsStacked =...
+          maskSubdivisionsStacked...
+          | maskSubdivisions(:,:,i) - bwmorph(maskSubdivisions(:,:,i), "remove");
+      end
+      output = maskSubdivisionsStacked;
+    end
+  end
+
+end
\ No newline at end of file