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

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

[Image Analyst]

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]

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.

"Pacman shape reference" was based on last question [link]. I changed the shape as follows:

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?

Answer 1

[Matt J]

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)

[Matt J]

Or you can use inpolygon,

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

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

mask=inpolygon(X,Y,xy(:,1),xy(:,2));

imshow(mask)

[waqas]

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]

But then the overall shape you are trying to achieve is just a combination of polygons and pacmen. Since you already know how to create each one separately, I do not see why you cannot just create an individual mask for each and or() them together.

[waqas]

I mean getting a mask is not the only goal. Idea is to get theta and radius of the grid points that would lie inside the mask region. So using a combination of pol2cart and cart2pol is something that should give the desired output variable but, I am unable to figure this out. The initial code that I have included in the main code of this question is doing something like this but the mask is not exactly what I am trying to get.

[Matt J]

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]

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]

Maybe you need

[y,x]=find(mask)

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

[waqas]

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]

Well, I know nothing about what is being computed across the mask, and so no grasp on why certain regions should be positive or negative. Can't you just take a sample point in the green region and check to see if its polar coordinates are being correctly computed? If the coordinates are wrong, then it is easier to start investigating why they are wrong.

[waqas]

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]

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]

Yeah this is precisely what I want to implement now. To add the translation vector to xt1 and yt1. But, I don’t know how to do that with grid data. Can you help with that?

[Matt J]

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]

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