How to find a constrained sum within a region

5 views (last 30 days)
Nate F on 30 Sep 2020
Commented: Image Analyst on 1 Oct 2020
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

Xingwang Yong on 30 Sep 2020
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 on 1 Oct 2020
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!

Image Analyst on 30 Sep 2020
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.
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.
Image Analyst on 1 Oct 2020
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.
% Get the area
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);