Image analysis - Bounding Box converted to Circular Mask

Hello, I have a greyscale image that I have located the centroids, binarised, area filtered and then drawn bounding boxes around

I create the bounding boxes by:

    s = regionprops(img1,'Centroid','Circularity','BoundingBox');
    cen = cat(1,s.Centroid);  circ=cat(1,s.Circularity);
    x1=cen(:,1); y1=cen(:,2);
    BB=cat(1,s.BoundingBox);

so BB has all the info in.

As I want to project this mask onto my raw image to extract the integrated intensity at each spot, I would rather circularise the boxes. This has been done for a single spot using:

     % For e.g the 25th spot  
     i=25;
     R=max(BB(i,3),BB(i,4));     % take the max dimension of the rectangle
     V=nsidedpoly(1000,'Center',[x1(i) y1(i)],'Radius',R*1).Vertices;  
     mask=poly2mask(V(:,1),V(:,2),sy,sx);
     ax4=nexttile; 
     myImagesc(app,ax4,mask);  title(ax4,'Keep');

This works beautifully for a single spot

So my question is how I get all of the spots on the binary mask , I tried the following but with no luck:

    n=numel(x1);        % This is the number of spots
    mask=[];
    for i=1:n     
       R=max(BB(i,3),BB(i,4));
       V=nsidedpoly(1000,'Center',[x1(i) y1(i)],'Radius',R*1).Vertices;  %Was 10
       mask=mask+poly2mask(V(:,1),V(:,2),sy,sx);
    end    

I dont know how to correctly perform the "addition or combining" of all the individual circular masks

    mask=mask+poly2mask(V(:,1),V(:,2),sy,sx);

Answer 1

[Image Analyst]

To display circles rather than bounding boxes, you can just ask regionprops for the Centroid and EquivDiameter. Then use viscircles. Untested code:

props = regionprops(mask, grayscaleImage, 'Centroid', 'EquivDiameter', 'Area', 'MeanIntensity');
centers = vertcat(props.Centroid);
radii = [props.EquivDiameter] / 2;
% Get integrated intensity over each blob
integratedIntensity = [props.MeanIntensity] .* [props.Area];
% Display circles around centroids.
hold on;
viscircles(centers, radii, 'Color', 'r', 'LineWidth', 2);

% Find the blob with the highest integrated intensity
maxIntIntensity = max(integratedIntensity);
% Find which blob(s) have that
indexes = find(integratedIntensity == maxIntIntensity);
% Put a crosshairs over them
for k = 1 : numel(indexes)
    thisIndex = indexes(k);
    plot(centers(thisIndex, 1), centers(thisIndex, 2), 'r+', 'MarkerSize' = 40, 'LineWidth', 3);
end

[Jason]

Oh thats awesome - thanks I.A!

Answer 2

[Jason]

This seems to work:

            result=[];
            for i=1:n     
               R=max(BB(i,3),BB(i,4));
               V=nsidedpoly(1000,'Center',[x1(i) y1(i)],'Radius',R*1).Vertices;  %Was 10
               mask=poly2mask(V(:,1),V(:,2),sy,sx);
               if i==1
                  result=mask
               else
                  result(mask~=0)=true; 
               end
            end