**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Taking 1d cuts of data along polar contours

7 views (last 30 days)

Show older comments

Anshuman Pal
on 11 Nov 2019

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

### Accepted Answer

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

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 Comments

Anshuman Pal
on 11 Nov 2019

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
on 11 Nov 2019

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
on 11 Nov 2019

Edited: Adam Danz
on 11 Nov 2019

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
on 11 Nov 2019

Edited: Adam Danz
on 11 Nov 2019

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
on 12 Nov 2019

Perfect! Thank you very much.

Anshuman Pal
on 30 Nov 2019

@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
on 30 Nov 2019

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.

Anshuman Pal
on 30 Nov 2019

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

m = matfile('data_t3.mat');

data = m.data_t3;

Adam Danz
on 30 Nov 2019

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.

Anshuman Pal
on 4 Apr 2020

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
on 4 Apr 2020

Edited: Adam Danz
on 10 Apr 2020

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
on 4 Apr 2020

Anshuman Pal
on 4 Apr 2020

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

Anshuman Pal
on 9 Apr 2020

Edited: Anshuman Pal
on 9 Apr 2020

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

Anshuman Pal
on 9 Apr 2020

Edited: Anshuman Pal
on 9 Apr 2020

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
on 9 Apr 2020

Edited: Adam Danz
on 9 Apr 2020

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
on 10 Apr 2020

Adam Danz
on 10 Apr 2020

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
on 10 Apr 2020

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
on 10 Apr 2020

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

### More Answers (0)

### See Also

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)