How can I overlay a rotated grid onto a map and calculate the area-weighted properties

I have a map of some metric, shown as different colors in the figure above. Onto that map I want to overlay a grid which is rotated at some angle and may or may not have the same pixel size as the underlying background map. I have tried to visualize the grid by the two red rectangles which are tilted by an angle of 45 degree in the picture above. I want to know the area-weighted values of the underlying background map contained by each of the two rectangles or actually the overlying grid which consists of multiple of these rectangles.

The approach I have taken so far is to represent the background map as a vector with the corners of rectangles and to search in a nested loop for the intersections of the red target rectangles with these using the polybool function like this:

%loop over two target rectangles
for j=1:2
  %loop over background grid
  for i=1:n*n
      %X and Y are vectors representing the corners of the background grid, 
      %x and y vectors representing the corners two target rectangles (separated by NaNs)
      [xi, yi] = polybool('intersection', X(i,:), Y(i,:), x(j), x(j);
      area=polyarea(xi,yi);%partial intersection area
      weighted_value(j)=weighted_value(j) + grid_value(i)*area;%area-weighted values
      cum_area(j)=cum_area(j)+area;%calculate total area
  end
end
weighted_value=weighted_value./cum_area;%area-weighted average

Given that both my underlying map and overlaid grid may have thousands of pixels, this approach is however overly slow and I have thus tried to come up with a solution which avoids the two nested loops:

%X and Y here are cells, x and y vectors, but otherwise have the same meaning as above
[xi, yi] = polybool('intersection',X,Y,x,y);
xi=cell2mat(xi);%need to convert back to vector to be used by polyarea
yi=cell2mat(yi);
area=polyarea(xi.',yi.');%calculate area

While the above code works and is faster compared to the loop by a factor of 10, I am unable to retrieve the partial areas and the corresponding values of the underlying map as I did with the loop.

Does anybody have a suggestion how I could retrieve this information with a numerically efficient solution or maybe somebody can suggest a totally different approach to my problem?

Thanks, Georg