5 views (last 30 days)

Show older comments

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.

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.

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.

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.

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);

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!