SOS! PLEASE SAVE MY CODE:

In my code I have the following issues: 1. the count measurement is not accurately outputting the count of blobs in the zones of the image. 2. When I normalize the count the count becomes decimal values, how can I fix this? 3. I want to make sure that the code is not measuring any of the large black areas in the image. Can someone please help me fix this, thanks!

[Walter Roberson]

We are not going to read through hundreds of lines code mentally modeling everything that could possibly go wrong. We need some of your images, and you have to describe the difference between what you get and what you want for each one.

Answer 1

[Kevin Holly]

Instead of trying to figure out how you calculated the areas. I used the following approach to find the masks:

fig = figure;
fontSize = 12;

% directoryInfo = dir('*.tif');
allFileNames = {'C1_L1_WDa.jpg'};
numfiles = length(allFileNames);

% Define the pixel size in square microns
pixelSize = 1;
numZones = 12;

for m = 1:numfiles
    fprintf('Processing file #%d of %d : "%s".\n', m, numfiles, allFileNames{m});
    MyRGBImage = fullfile(pwd, allFileNames{m});
    % Read in image.
    imageData = imread(MyRGBImage);

    [rows, columns, numberOfColorChannels] = size(imageData);
    if numberOfColorChannels > 1
        grayImage = rgb2gray(imageData); % Convert to color.
    end

    % Apply background subtraction.
    grayImage = imtophat(grayImage, strel('disk', 100));
    %the lower this number the less lipid droplets
    % Apply contrast correction.
    grayImage = imadjust(grayImage);

    % Apply Gaussian filtering.
    grayImage = imgaussfilt(grayImage, 3);

    % Apply thresholding.
    thresholdValue = graythresh(grayImage);
    binaryImage = imbinarize(grayImage, thresholdValue);

    % Remove black elongated shapes.
    binaryImage = bwareaopen(binaryImage, 200);
    binaryImage = imclearborder(binaryImage);
    binaryImage = imfill(binaryImage, 'holes');

    % mask
    % Threshold the image to create a binary mask
    gray_img = rgb2gray(imageData); % Convert to color.
    threshold = 5;
    mask = gray_img > threshold;
    % Remove small objects
    min_size = 1000;
    mask = bwareaopen(mask, min_size);
    % Fill holes
    mask = imfill(mask, 'holes');
    % Erode and dilate the mask to remove noise and smooth edges
    se = strel('disk', 5);
    mask = imerode(mask, se);
    mask = imdilate(mask, se);
    
    close; 
    figure;
    % ZZ: Right now for each single image, a separate window (figure) 
    % will be created to plot the related info, and the former figure will
    % be closed.
    subplot(3, 4, 4);
    imshow(mask);
    title('Functional Mask, thres = 5');

    % Display the organelle segmented image.
    subplot(3, 4, 2);
    imshow(binaryImage, []);
    axis on;
    title('Organelle Segmented Image', 'FontSize', fontSize);

    % Display the original image.
    subplot(3, 4, 1);
    imshow(imageData, []);
    axis on;
    title('Original Image', 'FontSize', fontSize);
    % Set up figure properties:
    set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
    sgtitle(sprintf('Analysis of Organelles in Image "%s"', allFileNames{m}), 'FontSize', 28);

% %     % Locate the center
% %     promptMessage = sprintf('Click at the center of the periportal region');
% %     titleBarCaption = 'Continue?';
% %     buttonText = questdlg(promptMessage, titleBarCaption, 'OK', 'Quit', 'OK');
% %     if strcmpi(buttonText, 'Quit')
% %         return;
% %     end
% %     [x,y] = ginput(1);
    x = 4592.86363636364;
    y = 705.954545454545;
    
    % Find out what the max distance will be by computing the distance to each corner.
    distanceToUL = sqrt((1-y)^2 + (1-x)^2);
    distanceToUR = sqrt((1-y)^2 + (columns-x)^2);
    distanceToLL = sqrt((rows-y)^2 + (1-x)^2);
    distanceToLR= sqrt((rows-y)^2 + (columns-x)^2);
    maxDistance = ceil(max([distanceToUL, distanceToUR, distanceToLL, distanceToLR]));
    
    % Calculate the radius of each zone
    radius = linspace(0, maxDistance, numZones+1);
    
    subplot(3, 4, 3);
    imshow(mask, []);
    axis on;
    title('Bin drawing', 'FontSize', fontSize);
    hold on;
    line([x, x], [1, rows], 'Color', 'r', 'LineWidth', 2);
    line([1, columns], [y, y], 'Color', 'r', 'LineWidth', 2);

I added this section below

    % Preallocate
    zonemask = zeros(size(mask,1),size(mask,2),numZones+1);

    fig2 = figure;
    imshow(mask)
    figure(fig2)
    c(1) = drawcircle("Center",[x, y],"Radius",radius(1),"LineWidth",1,"Color","b","Visible","off");
    % Draw the circles for each zone
    for ii = 2:numZones+1
        figure(fig2)
        c(ii) = drawcircle("Center",[x, y],"Radius",radius(ii),"LineWidth",1,"Color","b","Visible","off");
        zonemask(:,:,ii) = createMask(c(ii))-createMask(c(ii-1));
        viscircles([x, y], radius(ii), 'LineStyle', '--', 'LineWidth', 1);

        total_intensity(ii) = sum(sum(gray_img.*uint8(mask).*uint8(zonemask(:,:,ii)))); % or replace binaryImage with mask
        total_zone_pixels(ii) = sum(sum(uint8(mask).*uint8(zonemask(:,:,ii))));
    end

    mean_intensity_zone = total_intensity./total_zone_pixels
end

Processing file #1 of 1 : "C1_L1_WDa.jpg".

mean_intensity_zone = 1×13

NaN 20.2110 18.5171 18.3373 16.1967 20.3758 21.0019 24.8756 23.1034 23.8200 22.3985 23.9102 16.3538

Verify Zones

figure
tiledlayout(5,3)
nexttile
imshow(zonemask(:,:,1))
nexttile
imshow(double(mask).*zonemask(:,:,1))
nexttile
imshow(gray_img.*uint8(mask).*uint8(zonemask(:,:,1)))

for ii = 2:5
    nexttile
    imshow(zonemask(:,:,ii))
    nexttile
    imshow(double(mask).*zonemask(:,:,ii))
    nexttile
    imshow(gray_img.*uint8(mask).*uint8(zonemask(:,:,ii)))
end

Below I made a colored version with the original image.

figure
mask2(:,:,1) = uint8(mask).*uint8(zonemask(:,:,ii));
mask2(:,:,2) = mask2(:,:,1);
mask2(:,:,3) = mask2(:,:,1);
imshow(imageData.*mask2)

[Kevin Holly]

As for counting the blobs, did you want to count the blobs within the zone based off of the image "mask" or "binaryImage"?

Once decided, multiply the zonemask by it and use regionprops.

rp = regionprops(mask.*zonemask) %or binaryImage instead of mask

The total amount of blobs within the zone would be

blobsinzone = length(rp)

Note, if you wanted the latter, you could replace "mask" in the above code with "binaryImage".

[Chanille]

Thank you so much for this I’m going to try to run it now.

[Chanille]

@Kevin Holly

Actually, I just have a question. What would be the difference in MATLAB in counting the zones based on the mask vs the binary image? I know they are different but the masks have more information than the binary. Would it just depend on what it is I am trying to measure in the image that defines which image to multiply by? Thank you.

[Kevin Holly]

Yes, if you were interested in finding the sum intensity values of all the blobs within the zone, I would use the binary image. Notice how I calculated the mean. I summed the grayscaled image after multiplying the masks and then divided by the total pixels of the blobs in the zone which would be the sum of matrix generated by the dot product of the binary mask and the zone mask.

I did not use the mean function as that would include the pixel values of zero (i.e. black).

[Chanille]

Oh I see now. I will double check this with the formula I am using to calculate the sum intensity of the pixels in my code. (Which by sum intensity I presume you mean the number of green blobs in the original image.)

[Kevin Holly]

Yes, the green blobs.

[Image Analyst]

I'm worried we have an XY Problem here. I'm not sure why you're segmenting the green blobs. If it's a fluorescent image and you want to know how much activity/fluorescence is in each annular zone you should just take the mean of the whole zone.

totalFluorescenceInThisZone = mean(gray_img(zoneMask));

If you're trying to threshold to avoid the black, don't worry about it. Its sum will be close to zero so it won't matter if it's included. But you can't just consider the sum since your annular masks all have different areas. You can't compare a zone with more area to one with less area so that's why you take the mean instead of the sum - to normalize out the different areas.

[Chanille]

@Image Analyst

What do you mean by XY problem sorry? Also here is how i am calculating the blobs. Is this correct? @Kevin Holly

 numObjects = length(stats);
        if numObjects > 0
            areas = [stats.Area];
            profileCounts(k) = length(stats)%-sum(index,"all");%# of objects should be numObjects but not 1
            totalArea(k) = nansum(areas)*pixelSize^2; %mc:this is supposed to be the total area of each zone
            avgSize(k) = nanmean(areas)*pixelSize^2;
            zoneArea(k) = areaOfEachZone(k); %mc:this is supposed to be the total area of each blob in the zone being measured
            circularities = [stats.Circularity];%circularities(index) = nan;
            avgCircularity(k) = nanmean(circularities);
            ferets = [stats.MajorAxisLength];% ferets(index)=nan;
            avgFeret(k) = nanmean(ferets)*pixelSize;
            minFerets = [stats.MinorAxisLength];% minFerets(index)=nan;
            avgMinFeret(k) = nanmean(minFerets)*pixelSize;
        end

[Walter Roberson]

@Chanille

An "XY Problem" is when a person asks to do something (X) because they think it will help them reach their goals, but their goals turn out to be something else (Y) that turns out not to need to use X after-all.

Suppose for example you were asking how to get epoxy to stick to anodized aluminum, but it turned out that your real question was about how to fasten down the piece of anodized aluminum. You would (in this example) have decided that you needed to "epoxy" to solve your problem, but "epoxy" turned out to be a distraction, and if you had stated from the beginning what you needed to fasten to what, then perhaps someone could have suggested bolts or welding or crazy glue.

[Chanille]

@Walter Roberson

Thanks. My question was answered!