How can I vectorize this code to decrease computation time?

I have created a script that takes two instances of a DICOM structure (e.g. the same tumor at two different points in time) and rotates one instance of the structure about the other instance's center of mass, for the three rotation planes, in an effort to find the x,y, and z axis rotation that leads to the maximum overlap between the two structures.

At this stage the user defines the max rotational search limit and the step size. The max rotational search limit includes negative and positive degrees of rotation. So, for a rotation of 5 degrees with 1 degree step size, as defined by user input there are 1331 possible combinations of rotations:

x = -5 -4 -3 -2 -1 0 1 2 3 4 5

y = -5 -4 -3 -2 -1 0 1 2 3 4 5

z = -5 -4 -3 -2 -1 0 1 2 3 4 5

The problem is, as the user selects higher max rotation limits my code slows severely. I anticipate that I could rewrite some of the for-loops using vectorization, however I am having significant difficulty and any help would be greatly appreciated. The code/pseudocode looks like this:

for x_theta = -MAX_ROT : STP_SZ : MAX_ROT 
  for y_theta = -MAX_ROT : STP_SZ : MAX_ROT  
      for z_theta = -MAX_ROT : STP_SZ : MAX_ROT 
R = makehgtform('xrotate', x_theta, 'yrotate', y_theta, 'zrotate', z_theta);
R = R(1:3, 1:3);

 for j = 1:1:length(array containing points of each slice)  
 s =  points on slice j  - center of mass point;                     % Subtract center of mass from points so we can rotate about the center of mass                                        
 rotatedPoints{:, j} = (R * s);                                      % Multiply points by rotation matrix
 rotatedPoints{:, j} = rotatedPoints{:, j} + center of mass point;   % Add center of mass back to points
 end

 MASK = false(dimensions of 3D image volume);

 for k = 1 : 1 : length(rotatedPoints)

      slice = interp1(SCNINFO.start(3):SCNINFO.width(3):...
                    SCNINFO.start(3) + (SCNINFO.dimensions(3) - 1) ...
                    * SCNINFO.width(3), 1:SCNINFO.dimensions(3), ...
                    -rotatedPoints{1, k}(1, 3), interpolator, 0);

           if slice ~= 0

                    mask = poly2mask((Rotatedpoints_temp{1, k}(:,2) + (SCNINFO.width(2) * ...
                                      (SCNINFO.dimensions(2) - 1) + ...
                                        SCNINFO.start(2)) + SCNINFO.width(2)/2) / ...
                                         SCNINFO.width(2) + 1, ...
                                          (Rotatedpoints_temp{1, k}(:,1) + SCNINFO.width(1)/2 - ...
                                           SCNINFO.start(1)) / SCNINFO.width(1) + 1, ...
                                            SCNINFO.dimensions(1), ...
                                             SCNINFO.dimensions(2));
                                     end

                    % If the new mask will overlap an existing value, 
                    % subtract
                    if max(max(mask + MASK(:,:,slice))) == 2
                        MASK(:,:,slice) = ...
                            MASK(:,:,slice) - mask;

                    % Otherwise, add it to the mask
                    else
                        MASK(:,:,slice) = ...
                            MASK(:,:,slice) + mask;
                    end       
            end

overlapImage = reference structure mask & MASK;
numOverlapPixels = nnz(overlapImage);
rotatedResults = 2 * numOverlapPixels / (volume reference structure + volume structure);     
      end
  end
 end

After running with the MATLAB profiler most time is consumed during the poly2mask. This is because I am making a binary mask from every slice of points for every rotation combination!