**You are now following this question**

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

# Getting meshgrid along different angles for data extraction from a matrix/image?

4 views (last 30 days)

Show older comments

Hi,

I have three point locations of a matrix which I want to use to get inclined meshgrid about the lines joining these points. I know how to individually perform these but I am unable to figure out how to rotate both of them together using pol2cart function. A similar question where pacman shape is obtained using the same technique. However, in the current implementation, I do not want a pacman shape and the rounded region that we see in the second plot should be at the end of the second line.

How to resolve this?

x0 = 500; y0 = 1;

x1 = 450; y1 = 300;

x2 = 600; y2 = 700;

m = 1000; n=1500;

subplot(2,1,1)

plot([y0,y1,y2],[x0,x1,x2]); axis ij equal; xlim([1,n]);ylim([1,m]);

grid on;

title('This is how the lines should look in the image or matrix')

[yi,xi] = meshgrid(y1-(1:n),x1-(1:m));

ang = atan((y1-y0)/(x1-x0));

ang2 = atan((y2-y1)/(x2-x1));

[the,ro] = cart2pol(yi,xi);

tt = the(:,1:x1);rt = ro(:,1:x1);

[xt1,yt1] = pol2cart(tt+ang,rt);

tt2 = the(:,(x1+1):end);rt2 = ro(:,(x1+1):end);

[xt2,yt2] = pol2cart(tt2+ang2,rt2);

x = [xt1,xt2];

y = [yt1,yt2];

intR = 10; extR = 100;

subplot(2,1,2)

crack2 = x<=intR & x>=-intR & y<=sqrt(intR^2-x.^2);

notcrack= x<=(intR+extR) & x>=-(intR+extR) & y<=sqrt((intR+extR)^2-(x).^2);

region= notcrack & ~crack2;

imshow(region);h = gca; h.Visible = 'on';

title('this is the maximum external region')

##### 2 Comments

Image Analyst
on 22 Jul 2020

I don't understand. What is an "inclined meshgrid"?

And do you have some mathematical/analytical formula for the mask and you just want to create a mask image(s) for some region(s) like a mask for crack and a mask for non-crack?

Then you say "I do not want a pacman shape". Well, I do not see such a shape in the images you posted, so I'm not sure why you said that.

Is your "maximal external region" image what you want, or not? It has no pacman shape. If it's not what you want, then what do you want?

Then you say "the rounded region that we see in the second plot should be at the end of the second line." Well, I do not see a second plot. Just one plot and a binary image. Do you want the U-shaped region in the left of the image to somehow be on your graph of the two line segments in the plot? I don't understand what that would look like so please show us.

waqas
on 23 Jul 2020

Edited: waqas
on 23 Jul 2020

What I actually want to find out is the value of a function similar to:

The values are defined by where you take the crack tip and overall angle of the crack in the image. Value of n varies from -1 to 2. For a straight crack, I can inclined the mesh using pol2cart as was answered in last question I asked (link to open the question). But if the crack is changing directions e.g., in the first figure of current question where the crack tip is now at (600,700), then implementation of pol2cart should change. Thats why I used "inclined meshgrid".

I do not have a formula for mask but the idea is that the mask would consist of some region outside the crack too. Howoever, crack itself will be removed from the mask as the formula I mentioned above is not suppose to be applied on the crack itself.

As you can see that the round region is after the crack tip only.

"Maximum external region" in would be R_ext now and I want to make a mask that would take care of the change of direction of the crack, which is shown in the plot attached with the question.

Hopefully I have explained all the points clearly. Now, I want to select region around two line (in the plot) to be equal to R_ext and then at the end of the crack tip i.e., (600,700), I want to make a rounded shape just like the figure posted above. Eventually I am going to use this mask to get the pixel values of grayscale images that I have and then analyze them.

In my previous question, you mentioned that we can use image processing techniques to get the mask too. Can you share the details?

### Accepted Answer

Matt J
on 22 Jul 2020

If you're just trying to generate shaded polygons, it would be simpler just to use impoly, e.g.,

x0 = 1; y0=1;

x1 = 200; y1=500;

x2 = 500; y2=100;

m = 800; n=600;

xy=[x0,y0;x1 y1;x2 y2];

xy=[xy;flipud(xy+[0,100])];

H=impoly(ancestor(imshow(false(m,n)),'axes'),xy);

imshow(H.createMask)

##### 14 Comments

waqas
on 23 Jul 2020

Now, I am getting the mask from the following code:

n = 2064;m = 1088;

xv = [1 300 900 900 300 1 NaN 1 300 700 700 300 1];

yv = [710 660 810 400 250 300 NaN 500 450 600 610 460 510];

[X,Y]=meshgrid(1:n,1:m);

mask=inpolygon(X,Y,xv,yv);

imshow(mask)

But, as I replied to Image Analyst's comment, my end goal is to get radius and angle (with respect to crack tip) values of the points which are in the mask region to implement my fomula. Crack tip in this figure is the point where the black line inside white region is ending. I have provided more details of how it looks based on modification of code from previous question which had just a straight crack.

Matt J
on 23 Jul 2020

waqas
on 23 Jul 2020

Matt J
on 23 Jul 2020

Edited: Matt J
on 23 Jul 2020

Idea is to get theta and radius of the grid points that would lie inside the mask region.

Well, once you have the mask, you can use find() to get the Cartesian coordinates. Then, you could apply cart2pol directly.

[x,y]=find(mask);

[theta,rho] = cart2pol(x-x0,y-y0)

waqas
on 23 Jul 2020

Edited: waqas
on 23 Jul 2020

I implemented that but I think there is a problem with current approach. In the attached figure, I think now the yellow line is the taken as axis line around which the angles and radii are calculated. This would make the region that is highlighted as green wrong. Current implementation would mean that crack is starting from (310,1) and ends around (600,540). Eventually when we use the formula that I have posted in one of the comments above, this region would have wrong results. If not, then may be I am making a mistake in the implementation so I am sharing the code here too.

n = 2064;m = 1088;

xv = [1 300 900 900 300 1 1 NaN 1 300 700 700 300 1 1];

yv = [710 660 810 400 250 300 710 NaN 500 450 600 610 460 510 500];

[X,Y]=meshgrid(1:n,1:m);

mask=inpolygon(X,Y,xv,yv);

[y x] = find(mask);

[theta rho] = cart2pol(y-605,x-700);

imshow(mask)

b = NaN(size(X));

b(mask) = theta;

subplot(2,1,1)

imshow(b)

imshow(b,[]);h = gca; h.Visible = 'on';title('theta')

c = NaN(size(X));

c(mask) = rho;

subplot(2,1,2);imshow(c,[])

h = gca; h.Visible = 'on';title('rho')

Matt J
on 23 Jul 2020

Maybe you need

[y,x]=find(mask)

You need to be sure you are using consistent coordinate conventions throughout.

waqas
on 23 Jul 2020

But that does not resolve the problem about the green zone.

May be following figure would clear my point. In the figure, the crack is not changing direction i.e., slope is almost the same. So, I can apply the solution that you gave in previous question and get the desired solution. However, the crack in this question is not along the same line and if I apply the above mentioned mask then in the green region I will get dark blue (negative displacements) while the displacement should be positive i.e., at the magnitude of 15 according to the colorbar.

Matt J
on 23 Jul 2020

Edited: Matt J
on 23 Jul 2020

waqas
on 26 Jul 2020

Hi,

I have done some modifications to the code that I originally posted. I think with just a little bit of help, my problem would be resolved and I will be able to close the question too. I have added description the figure as well as the code to help you or others understand the steps that I am doing or logic behind them.

- What I want to do right now is to shift the red line for the first part of the crack to the first part of the blue line that I plotted in subplot(2,1,1). I could not figure out how to shift the grid. Rotation part is clear to me. We just need to shift the xt1 and yt1 accordingly.
- Hopefully after that we will see the mask for first half the crack as well.

Looking forward to you comments. [Figure and Code below]

x0 = 500; y0 = 1;

x1 = 400; y1 = 500;

x2 = 600; y2 = 1000;

m = 1088; n=2060;

subplot(2,1,1)

plot([y0,y1,y2],[x0,x1,x2]); axis ij equal; xlim([1,n]);ylim([1,m]);

grid on;

title('This is how the crack should look in the image or matrix')

[yi,xi] = meshgrid(y2-(1:n),x2-(1:m));

ang = atan((y1-y0)/(x1-x0)); %%The rotation that will be added to first part of the crack

ang2 = atan((y2-y1)/(x2-x1));%%The rotation that will be added to second part of the crack

[the,ro] = cart2pol(yi,xi);

tt = the(:,1:x1);rt = ro(:,1:x1); %Selecting region till first crack tip

[xt1,yt1] = pol2cart(tt+ang,rt); % adding the "ang" for rotating the initial grid to meet the crack line

tt2 = the(:,(x1+1):end);rt2 = ro(:,(x1+1):end); %Selecting region for second part of crack

[xt2,yt2] = pol2cart(tt2+ang2,rt2); %Adding "ang2" to get rotation along brown line in plot 1

x = [xt1,xt2]; %Combine x coordinates of part 1 of the crack with second part crack

y = [yt1,yt2];

intR = 10; extR = 200; %intR is the region from the inside, while extR is the maximum external region

subplot(2,1,2)

crack2 = x<=intR & x>=-intR & y>=-sqrt(intR^2-x.^2); % to remove the internal region

notcrack= x<=(intR+extR) & x>=-(intR+extR) ; %For external region

notcrack2 = y>=-sqrt((intR+extR)^2-(x).^2);

region= notcrack & ~crack2 & notcrack2; % Final region

%%TODO: How to get both parts of the crack after shifting

imshow(region);h = gca; h.Visible = 'on';

title('Mask')

Matt J
on 26 Jul 2020

I don't understand too much of the new problem (you should probably post a more general explanation in a new thread), but if the question is how to find the translation vector that will put the red line tip to tip with the leftmost blue segment, then it is just the vector difference between the two tips:

translation_vector=[500,400]-[0,800]

waqas
on 26 Jul 2020

Matt J
on 26 Jul 2020

I may be missing something, but it seems to me it would just be,

xt1 = xt1 + translation_vector(1);

yt1 = yt1 + translation_vector(2);

waqas
on 26 Jul 2020

Thank you the discussion was really helpful. I am posting the complete code with all the corrections in case some needs help. I will accept your answer to close the thread too.

This is how the final mask looks like:

There is still a little inconsistancy right where the lines meet but I can live with that.

x0 = 500; y0 = 1;

x1 = 400; y1 = 500;

x2 = 600; y2 = 1000;

m = 1088; n=2060;

subplot(2,1,1)

plot([y0,y1,y2],[x0,x1,x2]); axis ij equal; xlim([1,n]);ylim([1,m]);

grid on;

title('This is how the crack should look in the image or matrix')

[yi,xi] = meshgrid(y2-(1:n),x2-(1:m));

% ang = atan((y1-y0)/(x1-x0));

% ang2 = atan((y2-y1)/(x2-x1));

ang = atan((x1-x0)/(y1-y0));

ang2 = atan((x2-x1)/(y2-y1));

[the,ro] = cart2pol(yi,xi);

tt = the(:,1:y1);rt = ro(:,1:y1);

[yt1,xt11] = pol2cart(tt-ang,rt);

translation = (x0-x2) + ((y2-y0)*tan(ang));

xt1 = xt11 + translation;

tt2 = the(:,(y1+1):end);rt2 = ro(:,(y1+1):end);

[yt2,xt2] = pol2cart(tt2-ang2,rt2);

x = [xt1,xt2];

y = [yt1,yt2];

intR = 10; extR = 200;

subplot(2,1,2)

crack2 = x<=intR & x>=-intR & y>=-sqrt(intR^2-x.^2);

notcrack= x<=(intR+extR) & x>=-(intR+extR) ;

notcrack2 = y>=-sqrt((intR+extR)^2-(x).^2);

region= notcrack & ~crack2 & notcrack2;

imshow(region);h = gca; h.Visible = 'on';

title('Mask')

### More Answers (0)

### See Also

### Categories

### 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 (한국어)