How to find a constrained sum within a region

Hi all! This might be a simple question and I do already have a method but I'm wondering if anyone has a better way to do this.

The problem is how to find the point of a 2D distribution quickly. Lets say you have a 2D Gaussian in an NxN matrix Z fromed from where w is the width, and X and Y are both matricies formed form meshgrid. The point of the Gaussian can also be thought of where 86% of the data lies within the function.

One way to find the point is to normalize the matrix Z to the sum of all the elements and then start summing all elements from the center of the distribution until the condition of 86% of the data is enclosed. The example code is below.

My question is - can this be done in a more efficient manner and if so can the resolution of the steps to enclose the point be finer than the grid size created via meshgrid?

S = 0;
Step_in = 0.001;
index = 1;
while S < 0.86                                              % Stay in Loop until S encloses 86% of the data
    % Data is the normalized distribution matrix
    Region = circ(X,Y, Step(index));                        % Circ is a function that creates a circular region of Diameter Step(index)
    S(index + 1) = sum(sum(Region.*Data))/sum(sum(Data));   % Find the sum of elements enclosed by the circular region
    Step(index+1) = Step(index) + Step_in;                  % continue making a larger circular region until the loop condition is met
    index = index + 1;
end

Answer 1

[Xingwang Yong]

You can do this in a binary search manner, instead of doing sum() stepwise.

Below is pseudo code

diameterMin = 0.001; % choose a diameter min that make S<0.86
diameterMax = 1; % choose a diameter max that make S>0.86
diameterMid = (diameterMin + diameterMax) / 2;
S = 0;
while abs(S - 0.86) > 0.001 
    S = calculate_S(diameterMid);
    if S > 0.86
        diameterMax = diameterMid;
    else
        diameterMin = diameterMid;
    end
    diameterMid = (diameterMin + diameterMax) / 2;
end


function S = calculate_S(diameter)
% ...
end

As for the resolution, I do not think you can have a finer resolution than gride size.

[Nate F]

I just tried this method! It seems to work a lot faster, it doesn't seem to always find the right point and it gets stuck in an infinite loop. I shall look into this to see if I can get it to work more consistently. Thanks for the idea!

Answer 2

[Image Analyst]

Not sure what you mean by "the 1/e^2 point" but try this and see if it's what you want.

numRows = 600; % Whatever you want.
[x, y] = meshgrid(1:numRows, 1:numRows);
w = 200;
z = exp(-(x(:) .^ 2 + y(:) .^ 2) / w^2);
zImage = reshape(z, numRows, numRows);
imshow(zImage, []);
% Sort z
sortedZ = sort(z, 'descend');
% Get the number of pixels that will be 1/e^2
threshold = 1 - exp(-2) % This is 86.46% of the pixels in the image.
% Find the value for that
index = round(numel(z) * threshold)
zThreshold = sortedZ(index)
% Threshold the image.
mask = zImage >= zThreshold;
boundaries = bwboundaries(mask);
b = boundaries{1};
xb = b(:, 2);
yb = b(:, 1);
% Plot the zone covering 1/e^2 of the area.
hold on;
plot(xb, yb, 'r-', 'LineWidth', 4);
caption = sprintf('Pixels contained in red outline = %.3f%% of the image', 100*threshold);
title(caption, 'FontSize', 20);

So I'm giving you the (x,y) coordinates of the boundary that contains 86% of the value, and the intensity (z value) where that occurs. If that's not what you want, then explain better.

[Thomas Castle]

Appologies for hijacking your answer but you seem to be a god in image processing, would it be possible to get your advice on this problem?

https://www.mathworks.com/matlabcentral/answers/602512-identification-of-human-upper-body-from-photo-and-overlaying-an-image

Thanks!

[Nate F]

Thanks so much for your reply!

Let me clarify a bit further. The area enclosed that I am looking for is where 86% of the intensity lies not 86% of the pixels. Does that make sense? Let me give a complete set of code that does what I was saying. Maybe this will make it more clear as to what I mean by 86% of the intensity values. The code I showed here does exactly what I want but it gets absurbly slow, so I was wondering if there is another way to go about this to speed it up.

clear;clc;close all;
N = 512; dx = 1e-3;
[X,Y] = meshgrid((-N/2:N/2-1)*dx);                      % Create Grid Centered on (0,0)
w0 = 0.005:0.005:0.1;                                   % Arbitrary shape standard deviation of a Guassian

for NN = 1:length(w0)
    Data = exp(-2*(X.^2 + Y.^2)/w0(NN)^2);              % Specific Gaussian Distribution
    
    index = 1; S = 0;                                   
    Step(1) = w0/2;                                     % Set the first step size of the circular mask
    Step_in = dx/4;                                     % Radius step size
    while S < 0.86
        Region = circ(X,Y, Step(index));                % Circular Region
        S(index + 1) = sum(sum(Region.*Data))/sum(sum(Data));   % Summing up all points to see how much of the intensity matrix lies within a region
        Step(index+1) = Step(index) + Step_in;                  % Increase radius size if loop condition isn't met
        index = index + 1;
    end
    Interp_pts = 0:dx/100:Step(end);                            % Create a finer grid size
    vq = interp1(Step,S,Interp_pts) - 0.86;                     % Interpolate the radius size matrix to find a zero crossing
    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);        % Code to find a singular zero crossing by detecting where the first sign change happens
    fwhm = Interp_pts(zci(vq));                                 % Finding the first point of zero crossing which becomes the w0 value
    T(NN) = fwhm;
  NN
end
figure(1);          % Plot to check results, they should be equal to w0
plot(w0,T/2,'-bo'); hold on;
plot(w0,w0,'-rx');
title('Checking the results of finding the 86% point')
xlabel('Input standard deviation (w0)')
ylabel('Calculated w0')

function C = circ(x,y,D)
r = sqrt(x.^2+y.^2);
C = double(r<D/2);
C(r==D/2) = 0.5;
end

[Image Analyst]

Not sure why you need to do all that complicated stuff when you can just create the array and threshold it to compute the area:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;

numRows = 600; % Whatever you want.
[x, y] = meshgrid(1:numRows, 1:numRows);
w = 200;
z = exp(-(x(:) .^ 2 + y(:) .^ 2) / w^2);
zImage = reshape(z, numRows, numRows);
imshow(zImage, []);
% Get the number of pixels that will be 1/e^2
zThreshold = max(z) * (1 - exp(-2)) % This is 86.46% of max intensity in the image.
% Threshold the image.
mask = zImage >= zThreshold;
% Get the area
area = sum(mask(:));
boundaries = bwboundaries(mask);
b = boundaries{1};
xb = b(:, 2);
yb = b(:, 1);
% Plot the area within zThreshold of the max intensity of the image.
hold on;
plot(xb, yb, 'r-', 'LineWidth', 4);
caption = sprintf('Red outline at %.3f of the max intensity of the image.\nArea = %d pixels.', zThreshold, area);
title(caption, 'FontSize', 12);