Clear Filters
Clear Filters

Taking 1d cuts of data along polar contours

8 views (last 30 days)
I have 2d data saved in a matrix. I want to take 1d cuts of this along certain contours. However, since the data has a highly polar nature, the x-y matrix index coordinates are not very useful since they would only give me cuts along x=const. and y=const. lines. Instead, I would like to consider the (r,theta) polar coordinates, and take cuts of the data along circles and radial lines. How do I do this?
For example, consider the case of a 2d Gaussian:
xc = 0;
yc = 0;
sigma = 1;
[X, Y] = meshgrid(-2:0.1:2,-2:0.1:2);
exponent = ((X-xc).^2 + (Y-yc).^2)./(2*sigma^2);
amplitude = 1 / (sigma * sqrt(2*pi));
% The above is very much different than Alan's "1./2*pi*sigma^2"
% which is the same as pi*Sigma^2 / 2.
z = amplitude .* exp(-exponent);
example_polar plot.png
So here, I would like to take 1d cuts of the data along the red circles (r=const.) and the blue lines (theta=const.). Please advise

Accepted Answer

Adam Danz
Adam Danz on 11 Nov 2019
Your data are quite coarse. You can see the unit squares that make up each coordinate. So, you won't be able to create a perfect circle from your data. You can, however, isolate the coordinates that are closest to the perimeter of a circle. Follow this code to reproduce the image below. The red x's mark coordinates that are closest to the circle perimeter.
%% Your code
xc = 0;
yc = 0;
sigma = 1;
[X, Y] = meshgrid(-2:0.1:2,-2:0.1:2);
exponent = ((X-xc).^2 + (Y-yc).^2)./(2*sigma^2);
amplitude = 1 / (sigma * sqrt(2*pi));
z = amplitude .* exp(-exponent);
%% My Addition
% Make aspect ratio equal
axis equal
% draw lines that pass through center of data
xCnt = 21; % X center; you can also do round(size(z,2)/2)
yCnt = 21; % Y center; you can also do round(size(z,1)/2)
% Define radius
r = 10;
% Draw circle over sampling area
hold on
th = linspace(0,2*pi,150); %100 is arbitrary but should be much finer res than your data
xCirc = r * cos(th) + xCnt;
yCirc = r * sin(th) + yCnt;
h = plot(xCirc, yCirc, 'k-');
% Your data are coarse but the circle perimeter is fine.
% Here we'll grab your data closest to the circle perimeter
[xDat, yDat] = meshgrid( 1:size(z,2),1:size(z,1));
pd = pdist2([xDat(:), yDat(:)], [xCirc(:), yCirc(:)]);
% Find the min dist for each circle sample (for each column)
[minDist, minDistIdx] = min(pd,[],1);
% Since the cirlce has a finer resolution than your data, remove consecutive duplicates
minDistIdx = minDistIdx([true,diff(minDistIdx)~=0]);
% Here are the x,y coordinates that are under the cirlce
xNearCirc = xDat(minDistIdx);
yNearCirc = yDat(minDistIdx);
% Plot those coordinates to confirm
plot(xNearCirc, yNearCirc, 'rx')
% Extract the z values of those coordinates
191111 155421-Figure 1.png
If this approach is suitable, you can also apply it to straight line.
Anshuman Pal
Anshuman Pal on 10 Apr 2020
Ah, I believe it's the diagonal. A final question, I believe -- what is the role of adjDist in the code? It's because I never understood that that I just accepted it without questioning.
Adam Danz
Adam Danz on 10 Apr 2020
Edited: Adam Danz on 10 Apr 2020
"But what I am looking for is 1d data... Do I need the diagonal?"
No. You need to convert the subscripts to a linear index.
linIdx = sub2ind(size(z), yNearCirc, xNearCirc)
zDat = z(linIdx);
and to confirm you did it correctly,
[xx,yy] = meshgrid(1:size(z,2), 1:size(z,1)); % replace with better variable names
plot(xx(linIdx), yy(linIdx), 'mo')
what is the role of adjDist in the code?
If you zoom into the coordinates of the black circle shown here, you'll see that they can fall anywhere within each square pixel. When you search for all pixels that are "close to" the circle's circumference in the line below, you need to +/- sqrt(2) from the radius in order to ensure that you're not eliminating pixels that contain the coordinates of the circumference. Why sqrt(2)? Your pixels are 1x1 and the length of the diagonal of a 1x1 square is sqrt(2).
idx = hypot(xDat(:)-xCnt, yDat(:)-yCnt) >= (r-adjDist) & ...
hypot(xDat(:)-xCnt, yDat(:)-yCnt) <= (r+adjDist);

Sign in to comment.

More Answers (0)


Community Treasure Hunt

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

Start Hunting!