Taking 1d cuts of data along polar contours

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

Answer 1

[Adam Danz]

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.

[Anshuman Pal]

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]

Also, what exactly is the second half of this chunk doing?

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

Should the curve always have finer resolution than the data for this to work? It's easy to do, of coures, but I just want to know what's going on.

[Adam Danz]

1. I see you've added the x and y inputs to imagesc(). That changes how the image is positioned within the plot. The xCnt and yCnt therefore needs new definitions.

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

2. That also changes the coordinates we're matching with the circle perimeter.

% Replace this line
[xDat, yDat] = meshgrid( 1:size(z,2),1:size(z,1));
% with these
xDat = X; % Defined earlier in your code
yDat = Y; % Defined earlier in your code

[Adam Danz]

Yes, the circle should ideally have finer sampling than your data. Here's why. Below is an image of the cirlce but I zoomed in and you can see each cirlce coordinates (I plotted with '.' marker) and I added grid lines so you can see your data units. For each black dot the code searches for the nearest square. As you can see, most squares contain more than 1 black dot. In order not to end up with a bunch of duplicate square selections, we eliminate all duplicates. That's what the 2nd line of code is doing in your question above.

You can imagine if the cirlce was sampled coarsly, it would skip over many squares. In fact, you can try that to see the results (fewer red X's). To decrease the sampling of the circle, reduce the integer here: th = linspace(0,2*pi,150);

[Anshuman Pal]

Perfect! Thank you very much.

[Adam Danz]

Glad I could help!

[Anshuman Pal]

@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?

[Anshuman Pal]

For reference, here is how it looks when I find the maximum of the original, coarse data:

The only lines I change are:

figure
imagesc([0 numCols], [0 numRows], Vq);
...
% find position of maximum
[max_num, max_idx] = max(data(:));
[y0, x0]=ind2sub(size(data),max_idx)
...
% Define radius 
r = 14;

Beyond this stage, I use the interpolated data for finding the cut on the circles.

[Adam Danz]

Any chance you could attach a mat file containing "data" and any other variables I need to run that code?

[Anshuman Pal]

Yes, the file is attached. All the extra code you should need is:

m = matfile('data_t3.mat');
data = m.data_t3;

[Adam Danz]

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

1. the center of the bivariate distribution is in the center of your plot
2. 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.

[Anshuman Pal]

Hi @Adam Danz,

I am following up on this question after a long time. Everything works fine, but I have problems scaling this method up to cases where the data matrix is very fine (say 200x200 elements). In such a case, pdist2() is forced to make a distance matrix which is huge, on the order of hundreds of GB. So, Matlab just throws up an error and refuses to proceed.

Can you think of a way of implementing your algorithm without creating such a huge distance matrix?

And again, I would like to thank you for your help with this.

[Adam Danz]

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]

PS, I think the solution above is sufficient but another approach to think about is using the contour() function to get circular contour lines and then using getContourLineCoordinates() from the file exchange to extract the controur coordinates.

[Anshuman Pal]

Your solution works perfectly! Thank you very much again :)

[Adam Danz]

Thanks for the feedback.

[Anshuman Pal]

Hi Adam, I have another follow-up question. Suppose the circlular contour I choose oversteps the boundaries of the data, like here:

The cyan crosses here represent [xDat, yDat] after running your code snippet above. Firstly, I see how the cyan can live outside the data limits, but I don't understand how the red circle, which represents `minDistIdx`, can live outside the data limits.

But that's not a problem per se. My real question is -- given this `minDistIdx`, how do I modify your code so that my final 1d data set along the contour has z(minDistIdx) inside the data limits, but zero (or NaN, or -1) outside? Thus, the the size of `minDistIdx` stays the same, which means I am still sampling the full circle

[Adam Danz]

Hint:

xDat(xDat<0) = 0; 
% or
xDat(xDat<0) = NaN; 
% or 
xDat(xDat<0) = [];

[Anshuman Pal]

Actually, there is another bug. I believe this is what was really bothering me above. With the new algorithm, minDistIdx refers to the indices in the new, reduced [xDat, yDat], not the original data. Thus, z(minDistIdx) gives me the wrong results. Do you know how to deal with this?

I have attached images of what is going in. ex2 is the correct data extraction along the red contour circle, while ex3 shows what I get when I use the cyan band.

[Adam Danz]

I don't know what the new algorithm is or how minDistIdx differs between the new algo and the version in my original answer.

What's the reduced [xDat, yDat] and how does it differ from z?

I don't know how the attached images relate to the imagesc(z) plot. Try again to describe the problem as simply as possible as if I don't know anything about the project.

[Anshuman Pal]

I apologise for my haste. I have tried to give a full explanation of my problem in the attached pdf. I hope it is sufficiently clear.

[Adam Danz]

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.

• getContourLineCoordinates
• extract-contour-data-from-contour-matrix-c

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]

Thank you very much. I am still going through this. Firstly, I apologise for a typo. By Vq,I meant the data z.

I have a question about what happens at the end. Suppose xNearCirc and yNearCirc are both of length 50. Then

zDat = z(yNearCirc, xNearCirc); 

is a matrix of size 50x50. But what I am looking for is 1d data, so it should just be an array of length 50. Which of the columns/rows in zDat is what I need? Do I need the diagonal?

[Anshuman Pal]

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]

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