Taking 1d cuts of data along polar contours

9 ビュー (過去 30 日間)
Anshuman Pal 2019 年 11 月 11 日

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);
imagesc(z) 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

サインインしてコメントする。

採用された回答

Adam Danz 2019 年 11 月 11 日
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.
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);
imagesc(z)
% 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)
xline(xCnt)
yline(yCnt)
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
z(minDistIdx) If this approach is suitable, you can also apply it to straight line.
25 件のコメント表示非表示 24 件の古いコメント
Adam Danz 2020 年 4 月 10 日
"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) & ...

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!