Taking 1d cuts of data along polar contours

Question

Anshuman Pal 2019-11-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/490486-taking-1d-cuts-of-data-along-polar-contours

编辑： Adam Danz 2020-4-10

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Adam Danz 2019-11-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/490486-taking-1d-cuts-of-data-along-polar-contours#answer_400985

在 MATLAB Online 中打开

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);
imagesc(z)
%% 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)
xline(xCnt)
yline(yCnt)
% 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
z(minDistIdx)

If this approach is suitable, you can also apply it to straight line.

25 个评论
显示 23更早的评论隐藏 23更早的评论

Anshuman Pal 2019-11-11

在 MATLAB Online 中打开

This is exactly what I needed. For the Gaussian above, I can make the data finer at will, so coarseness is not an issue.

However, when I do so, I get a bug with the plot itself. The imagesc() by itself runs fine, but when I add the rest of the code, I just get a white plot :

% finer mesh and rescaled imagesc
xc = 0;
yc = 0;
sigma = 1;
xmin = -2;
xmax = 2;
[X, Y] = meshgrid(xmin:0.1:xmax, xmin:0.1:xmax);
exponent = ((X-xc).^2 + (Y-yc).^2)./(2*sigma^2);
amplitude = 1 / (sigma * sqrt(2*pi));  
z = amplitude .* exp(-exponent);
imagesc([xmin xmax], [xmin xmax], z)
% Make aspect ratio equal
axis equal
% draw lines that pass through center of data
xCnt = round(size(z,2)/2);  % X center;  you can also do round(size(z,2)/2)
yCnt = round(size(z,1)/2);  % Y center;  you can also do round(size(z,1)/2)
xline(xCnt)
yline(yCnt)
% Define radius. Changed to have same scale as data.
r = 1; 
% 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')

How do I fix this?

Anshuman Pal 2019-11-30

在 MATLAB Online 中打开

@Adam Danz, I need some further help on the above question. I am using the code to find the maximum of some data and draw circles around it. If I use it directly with the original data, it works. However, since the original data is too coarse for me to find the exact maximum correctly, I want to do the whole process with the data interpolated over a finer mesh. However, this time, there is some bug in finding the squares closest to the circle. In the image below, the black circle marks where I want to measure it, while the red dots mark where it actually measures:

Here is the code:

[numRows, numCols] = size(data)
[X, Y] = meshgrid(1:numCols,1:numRows); % original grid
[Xq,Yq] = meshgrid(1:0.1:numCols,1:0.1:numRows); % fine grid, 10x finer here
Vq = griddata(X,Y,data,Xq,Yq,'cubic');
figure
imagesc(Vq);
% imagesc([0 numCols], [0 numRows], Vq);
title('Cubic Interpolation Over Finer Grid');
% find position of maximum
[max_num, max_idx] = max(Vq(:)); % use smoother interpolated data
[y0, x0]=ind2sub(size(Vq),max_idx)
% [max_num, max_idx] = max(data(:));
% [y0, x0]=ind2sub(size(data),max_idx)
% Make aspect ratio equal
axis equal
% draw lines that pass through center of data
xCnt = x0;  % X center;  
yCnt = y0;  % Y center;  
xline(xCnt);
yline(yCnt);
size(Xq) % NB: 1d contour needs to have much finer resolution than size(Xq)
% Define radius 
r = 140; % r=15 is too large
% Draw circle over sampling area
hold on
th = linspace(0,2*pi,2000); %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-');
% Use refined grid Xq, Yq, and interpolated data Vq from now on.
% Here we'll grab your data closest to the circle perimeter
xDat = Xq; 
yDat = Yq;
pd = pdist2([xDat(:), yDat(:)], [xCirc(:), yCirc(:)]); 
% Find the min dist for each circle sample (for each column)
[minDist, minDistIdx] = min(pd,[],1); 
% Since the circle 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 circle
xNearCirc = xDat(minDistIdx); 
yNearCirc = yDat(minDistIdx); 
size(minDistIdx)
% Plot those coordinates to confirm
plot(xNearCirc, yNearCirc, 'ro', 'MarkerSize',1)

Can you please help?

Adam Danz 2019-11-30

在 MATLAB Online 中打开

Hi Anshuman, the problem was very easy to identify by stepping through each line of your code in debug mode. It's as easy as placing a break point at the top of your code and pressing F10 to step through each line, line-by-line, until you notice a change in your plot.

There are two areas you need to fix.

Your plot is fine until you get to these two lines

xline(xCnt)
yline(yCnt)

then the plot goes all white. If you look at the values of xCnt and yCnt you'll notice that they are nowhere near the center of your data. When those lines are plotted, the rest of your data gets crushed and you can no longer see it. Those value are also causing problems when creating the xCirc and yCirc variables.

xCnt and yCnt should be the center of your bivariate distribution. The current calculations are fine as long as two assumptions are met

the center of the bivariate distribution is in the center of your plot
the image is scaled such that each block is 1 unit along the x and y axes.

Those assumptions were true with your original data but now that you're using xmin and xmax in imagesc(), assumption #2 is no longer true.

So, now we need to redefine xCnt and yCnt. If assumption #1 is always true, you can redefine those values as

xCnt = xmin + (xmax-xmin)/2;
yCnt = xmin + (xmax-xmin)/2;

The rescaling caused a second problem in this line below.

[xDat, yDat] = meshgrid( 1:size(z,2),1:size(z,1));

The purpose of this line and the next two lines is the find the distances between each coordinate in your black cirlce and the values in xDat and yDat which are coordinates of each little square in your bivariate plot. Now that your bivariate plot ranges from -2 to 2 we need to compute the position of each square differently. That line should be changed to

[xDat, yDat] = meshgrid(linspace(xmin,xmax,size(z,1)),linspace(xmin,xmax,size(z,2))); 

One final note, my answer wasn't meant to be a solution that would fit all data sets. It was meant to work with an example data set to show how you can apply it to your data. It's likely that you'll need to make additional changes which will require you to understand what the code is going. I recommend taking several minutes to check out the debugging link I provided above because that's the best way to step through your code and see what it's doing.

Adam Danz 2020-4-4

编辑：Adam Danz 2020-4-10

在 MATLAB Online 中打开

Hi Anshuman , glad to see your project is coming along.

Yes, pdist 2 can get overloaded very quickly. Here's a much more efficient solution. The main idea is to eliminate image pixes that are far away from the circumference before computing pdist2.

It uses the same variable names from my answer. See inline comments for details (some of which are important).

% Isolate data that are within a square around the circle.
% The +/- adjDist is to account for the circumference passing through the 
% corner of an image pixels.  **If your pixes aren't 1x1 this value should be adjusted!
adjDist = sqrt(2); 
% [xDat, yDat] = meshgrid( xCnt-r-adjDist : xCnt+r+adjDist, yCnt-r-adjDist : yCnt+r+adjDist); 
[xDat, yDat] = meshgrid( xCnt-r : xCnt+r, yCnt-r : yCnt+r); % Corrected 4-April-2020
% Compute the distance of each coordinate to the circle's center. 
% Then get rid of any coordinates that are far from circumference.  
idx = hypot(xDat(:)-xCnt, yDat(:)-yCnt) >= (r-adjDist) & ... 
    hypot(xDat(:)-xCnt, yDat(:)-yCnt) <= (r+adjDist);
% Eliminate unnecessary data
xDat(~idx) = []; 
yDat(~idx) = []; 
% Show which image pixels will be analyzed
hold on
plot(xDat, yDat, 'c+')
% The rest is unchanged (starting with this line)
pd = pdist2([xDat(:), yDat(:)], [xCirc(:), yCirc(:)]); 

The final result is this. Cyan + show pixels that were analyzed. Red x show pixels nearest to the circle.

Adam Danz 2020-4-10

在 MATLAB Online 中打开

The pdf was indeed helpful and the code was very familiar.

If you're extracting values over a contour you can use the one of the files below from the file exchange.

Before addressing the problem, I discovered a small error in the code. The error was easy to detect by zooming into the image and seeing that the cyan + and red x were not centered in the pixel. This is the kind of exploration you should do with your data, especially when you're using some else's code.

To fix that error, remove the adjDist terms from the meshgrid line below. After the correction, it should look like the line below. No other changes are needed to address that problem.

[xDat, yDat] = meshgrid( xCnt-r : xCnt+r, yCnt-r : yCnt+r);

"The problem is that, when I use the [xDat, yDat] returned by isolateDataNearContour() in dataCut_1d(), I get back a minDistIdx which no longer corresponds to the index positions on the original mesh! Thus, Vq(minDistIdx) gives back nonsensical results."

Vq isn't defined anywhere in the pdf nor anyone in this thread but I think I understand the problem. After you make the change the meshgrid described above, xNearCirc and yNearCirc will be integers that define the column and row indices of z.

Instead of using

z(minDistIdx)

replace that with

zDat = z(yNearCirc, xNearCirc);

To confirm that the z values make sense, you can add a colorbar to the plot, scale the colorbar to the range of values in z, and then add two colorbar ticks defining the range of zDat values.

zDat = z(yNearCirc, xNearCirc); 
cb = colorbar();
caxis([min(z(:)), max(z(:))])
cb.YTick = [min(zDat(:)), max(zDat(:))]

Anshuman Pal 2020-4-10

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 2020-4-10

编辑：Adam Danz 2020-4-10

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Taking 1d cuts of data along polar contours

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

25 个评论
显示 23更早的评论隐藏 23更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Taking 1d cuts of data along polar contours

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

25 个评论 显示 23更早的评论隐藏 23更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

25 个评论
显示 23更早的评论隐藏 23更早的评论