Effacer les filtres
Effacer les filtres

How to compute the ratio of the insection point in each grid on the boundary of a real structure

2 vues (au cours des 30 derniers jours)
In the Cartesian grid, the real structures are pixelized into the staircase approximation. To reduce the staircase error, the effective permittivity is computed, which is weighted by the fraction of the mesh step inside the material 1 or 2 (See the below figure). For instance, the structure is sphere. How I can compute the fraction ratio of four edges in each grid on the boundary of the sphere? Could you please give me some suggestions?
% build Cartesian grid
[Y1,X1]=meshgrid(y_1grid,x_1grid);
% build structure (sphere)
UCC= ((X1+dx/2).^2 + (Y1+dy/2).^2) <= (diam/2.0)^2;
eps_cc=eps(2)*(UCC==1)+eps(1)*(UCC==0);
% find the boundary
[Fx,Fy]=gradient(eps_cc);
boundary=(Fx~=0)|(Fy~=0);
% compute the ratio of intersection point
??
Big Thanks to you
  2 commentaires
Jiali
Jiali le 11 Sep 2019
Yes, I want to compute fraction for each pixel. That's why I try to subgrid the current one pixel into 10 smaller grids and determine whether each point is inside the circel or not.

Connectez-vous pour commenter.

Réponse acceptée

Bruno Luong
Bruno Luong le 9 Sep 2019
Modifié(e) : Bruno Luong le 9 Sep 2019
Assuming the boundary has one connexed piece
x=-3:12;
y=-3:10;
[X,Y]=meshgrid(x,y);
cx = 4;
cy = 3;
r = 5;
Z = (X-cx).^2 + (Y-cy).^2 - r^2;
C = contourc(x,y,Z,[0 0]);
Cx = C(1,2:end);
Cy = C(2,2:end);
[Cx; Cy] % fractional crossing is grouped here
close all
hold on
plot(X,Y,'-b');
plot(X',Y','-b');
plot(Cx,Cy,'-r');
plot(cx,cy,'ko');
axis equal
  2 commentaires
Jiali
Jiali le 16 Sep 2019
Your method is very smart. However, the key is to know the fraction of each pixel. Thus, based on the crossing points calculated by your method, I write the below code to compute the fraction for each pixel. Could you please give me some suggestions for the following code?
ratio_x=ones(length(x),length(y));
ratio_y=ratio_x;
for k=1:length(Cx)
for ii=1:length(x)
for jj=1:length(y)
if Cx(k)==x(ii) && (Cy(k)<=y(jj)+dy && Cy(k)>=y(jj))
if (x(ii)-cx)^2+(y(jj)-cy)^2-r^2<=0
ratio_y(ii,jj)=abs(y(jj)+dy-Cy(k))/dy;
else
ratio_y(ii,jj)=abs(y(jj)-Cy(k))/dy;
end
end
if Cy(k)==y(jj) && (Cx(k)>=x(ii) && Cx(k)<=x(ii)+dx)
if (x(ii)-cx)^2+(y(jj)-cy)^2-r^2<=0
ratio_x(ii,jj)=abs(x(ii)+dx-Cx(k))/dx;
else
ratio_x(ii,jj)=abs(x(ii)-Cx(k))/dx;
end
end
end
end
end
Bruno Luong
Bruno Luong le 16 Sep 2019
Modifié(e) : Bruno Luong le 16 Sep 2019
The pixels coordinates are
ii=interp1(x,1:length(x),Cx,'previous')
jj=interp1(y,1:length(y),Cy,'previous')
And the corresponding fractional are
ratio_x=interp1(x,1:length(x),Cx)-ii
ratio_y=interp1(y,1:length(y),Cy)-jj
Just loop on the list [ii; jj; ratio_x; ratio_y] and do whatever you need to do.
for pixel = [ii; jj; ratio_x; ratio_y]
i = pixel(1);
j = pixel(2);
rx = pixel(3);
ry = pixel(4):
% do something with i,j,rx,ry ...
end
Alternatively, to avoid using INTERP1 you can call contour with pixel coordinates
C = contourc(1:size(Z,2),1:size(Z,1),Z,[0 0]);
Cx = C(1,2:end);
Cy = C(2,2:end);
% these replace the INTERP1 stuffs
ii = floor(Cx);
jj = floor(Cy);
ratio_x = Cx-ii;
ratio_y = Cy-jj;
The process ii,jj,ratio_x,ratio_y

Connectez-vous pour commenter.

Plus de réponses (1)

darova
darova le 11 Sep 2019
'That's why I try to subgrid the current one pixel into 10 smaller grids and determine whether each point is inside the circel or not.'
THAT IS GREAT IDEA
See attached script

Catégories

En savoir plus sur MATLAB dans Help Center et File Exchange

Produits


Version

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by