Skip to content
Snippets Groups Projects
Commit 63febaaa authored by Tom Görtzen's avatar Tom Görtzen
Browse files

some minor fixes

the method my_triangle_fix.g contains all the main methods
eps is chosen as a global variable for now (to play around a bit)

example_tom_triangulation.g contains some examples
parent 0ab7dcdd
Branches
No related tags found
No related merge requests found
TwoTriangleIntersection:=function(coords_j,coords_k,eps)
local l;
local l,ccoords,x,y,normal,t_coords,t_normal,ints,different,Compared,orthog,c_coords,numb_int,not_on_edge,ints_points,o,p,res,c_normal,d1,diff,not_on_edges,possb_intersec,lambda,edge_param,vOedge,d2;
ccoords:=coords_j;
x := ccoords[2]-ccoords[1];
y := ccoords[3]-ccoords[1];
......
Read("my_triangle_fix.g");
Read("my_triangle_fix.g");
### Example for triangulation function (by Tom)
### Ikosahedron 2_1
### Icosahedron 2_1
# 2_1 Icosahedron coordinates
......
# returns true if the two lines lie on the same plane
OnSamePlane:=function(l1,l2)
local eps;
eps:=Float(10^(-8));
OnSamePlane:=function(l1,l2,eps)
if Norm2(Crossproduct(Crossproduct(l1[1]-l1[2],l1[1]-l2[2]),Crossproduct(l1[1]-l1[2],l1[1]-l2[1]))) < eps then
return true;
fi;
......@@ -10,28 +8,188 @@ end;
# Computes intersection of two Line Segments in 3D
# https://math.stackexchange.com/questions/270767/find-intersection-of-two-3d-lines
LineSegmentIntersection:=function(l1,l2)
local C,D,e,f,g,eps,sol,dotsol1,dotsol2;
eps:=Float(10^(-8));
LineSegmentIntersection:=function(l1,l2,eps)
local C,D,e,f,g,sol,dotsol1,dotsol2;
C:=l1[1];
D:=l2[1];
e:=l1[2]-l1[1];
f:=l2[2]-l2[1];
g:=D-C;
if OnSamePlane(l1,l2) and Norm2(Crossproduct(e,f)) > eps then
if Dot(Crossproduct(f,g),Crossproduct(f,e)) > eps then
if OnSamePlane(l1,l2,eps) and Norm2(Crossproduct(e,f)) > eps then
if Dot(Crossproduct(f,g),Crossproduct(f,e)) >= 0. then
sol := C + Norm2(Crossproduct(f,g))/Norm2(Crossproduct(e,f))*e;
else
sol := C - Norm2(Crossproduct(f,g))/Norm2(Crossproduct(e,f))*e;
fi;
dotsol1 := Dot(l1[1]-l1[2],sol-l1[2]);
dotsol2 := Dot(l2[1]-l2[2],sol-l2[2]);
if dotsol1 >= 0. and dotsol2 >= 0. and dotsol1 <= Dot(l1[1]-l1[2],l1[1]-l1[2]) and dotsol2 <= Dot(l2[1]-l2[2],l2[1]-l2[2]) then
if Norm2(l1[1]-sol)+Norm2(l1[2]-sol)<=Norm2(l1[1]-l1[2])+eps and Norm2(l2[1]-sol)+Norm2(l2[2]-sol)<=Norm2(l2[1]-l2[2])+eps then
if not NumericalPosition(l1,sol,eps)<>fail or not NumericalPosition(l2,sol,eps)<> fail then
return [true,sol];
fi;
fi;
fi;
return [false];
end;
#LineSegmentIntersectionColinear([ [ 1., 0., 0. ], [ 2., 0., 0 ] ],[ [ 1., 0., 0. ], [ 1.6, 0., 0] ],1.e-06);
#LineSegmentIntersectionColinear([ [ 1., 0., 0. ], [ 1.6, 0., 0 ] ],[ [ 1., 0., 0. ], [ 2, 0., 0 ] ],1.e-06);
#LineSegmentIntersectionColinear([ [ 1., 0., 0. ], [ 1.6, 0., 0 ] ],[ [ 1., 0., 0. ], [ 2, 0., 0 ] ],1.e-06);
#LineSegmentIntersectionColinear([ [ 1., 0., 0. ], [ 1.6, 0., 0 ] ],[ [ 1., 0., 0. ], [ 2, 0., 0 ] ],1.e-06);
#LineSegmentIntersectionColinear([ [ 1., 0., 0. ],[ 1.6, 0., 0 ]],[ [ 2., 0., 0 ] ,[ 1., 0., 0. ]],1.e-06);
#LineSegmentIntersectionColinear([ [ 1.6, 0., 1.2 ],[ 1., 0., 0. ] ],[ [ 1., 0., 0. ], [ 2., 0., 2. ] ],1.e-06);
# test for other cases:
choice:=[[ 1.6, 0., 0 ],[ 2., 0., 0 ]];
for i in [1..2] do
for j in [1..2] do
for k in [1..2] do
l1:=[];
l2:=[];
l1[i]:=[ 1., 0., 0. ];
l2[j]:=[ 1., 0., 0. ];
if BoundPositions(l1)[1]=1 then
l1[2]:=choice[k];
else
l1[1]:=choice[k];
fi;
if k=1 then
n:=2;
else
n:=1;
fi;
if BoundPositions(l2)[1]=1 then
l2[2]:=choice[n];
else
l2[1]:=choice[n];
fi;
#Print(LineSegmentIntersectionColinear(l1,l2,eps));
od;
od;
od;
LineSegmentIntersectionColinear:=function(l1,l2,eps)
local v1,v2,x1,x2,y1,y2,r1,r2,s1,s2,i,j;
l1:=1.*l1;
l2:=1.*l2;
v1:=l1[2]-l1[1];
v2:=l2[2]-l2[1];
if OnSamePlane(l1,l2,eps) and Norm2(Crossproduct(v1,v2)) < eps and PointsInOneLine(l1[1],l1[2],l2[1],eps) and PointsInOneLine(l1[1],l1[2],l2[2],eps) then
x1:=l1[1];
x2:=l2[1];
y1:=l1[2];
y2:=l2[2];
r1:=[];
r2:=[];
s1:=[];
s2:=[];
for i in [1..3] do
if Sqrt(v1[i]^2)>eps then
r1[i]:=(x2[i]-x1[i])/v1[i];
r2[i]:=(y2[i]-x1[i])/v1[i];
fi;
if Sqrt(v2[i]^2)>eps then
s1[i]:=(x1[i]-x2[i])/v2[i];
s2[i]:=(y1[i]-x2[i])/v2[i];
fi;
od;
# check if r1 has same non-zero entries
for i in BoundPositions(r1) do
for j in BoundPositions(r1) do
if r1[i]^2>eps and r1[j]^2>eps and (r1[i]-r1[j])^2>eps then
return [false];
fi;
od;
od;
for i in BoundPositions(r1) do
if r1[i]^2>eps then
r1:=r1[i];
break;
fi;
od;
if IsList(r1) then
r1:=0.;
fi;
for i in BoundPositions(r2) do
if r2[i]^2>eps then
r2:=r2[i];
break;
fi;
od;
if IsList(r2) then
r2:=0.;
fi;
for i in BoundPositions(s1) do
if s1[i]^2>eps then
s1:=s1[i];
break;
fi;
od;
if IsList(s1) then
s1:=0.;
fi;
for i in BoundPositions(s2) do
if s2[i]^2>eps then
s2:=s2[i];
break;
fi;
od;
if IsList(s2) then
s2:=0.;
fi;
if eps<= r1 and r1 <=1-eps then
if eps<=s1 and s1<=1-eps then
return [true,[[x1,x1+r1*v1],[x1+r1*v1,y2-s1*v2],[y2-s1*v2,x2]],1];
elif eps<=s2 and s2 <=1-eps then
return [true,[[x1,x1+r1*v1],[x1+r1*v1,x2+s2*v2],[x2+s2*v2,y2]],2];
fi;
fi;
if eps<=r2 and r2 <=1-eps then
if eps<=s1 and s1 <=1-eps then
return [true,[[y1,y1-r2*v1],[y1-r2*v1,y2-s1*v2],[y2-s1*v2,x2]],3];
elif eps<=s2 and s2<=1-eps then
return [true,[[y1,y1-r2*v1],[y1-r2*v1,x2+s2*v2],[x2+s2*v2,y2]],4];
fi;
fi;
# test if l1 contains l2 and a point of l2
if FlVEq(x1,x2,eps) then
#r1=0 and s1=0
if eps<=r2 and r2 <=1-eps then
return [true,[[x2,y2],[y2,y1]],5];
fi;
if eps<=s2 and s2<=1-eps then
return [true,[[x1,y1],[y1,y2]],6];
fi;
elif FlVEq(y1,x2,eps) then
#s2=0 and r1=1
if eps<=r2 and r2 <=1-eps then
return [true,[[x1,y2],[y2,x2]],7];
fi;
if eps<=s1 and s1<=1-eps then
return [true,[[y1,x1],[x1,y2]],8];
fi;
elif FlVEq(y2,x1,eps) then
#r2=0 and s1=1
if eps<=r1 and r1 <=1-eps then
return [true,[[y2,x2],[x2,y1]],9];
fi;
if eps<=s2 and s2<=1-eps then
return [true,[[x2,y1],[y1,x1]],10];
fi;
elif FlVEq(y2,y1,eps) then
#r2=1 and s2=1
if eps<=r1 and r1 <=1-eps then
return [true,[[x1,x2],[x2,y2]],11];
fi;
if eps<=s1 and s1<=1-eps then
return [true,[[x2,x1],[x1,y1]],12];
fi;
fi;
fi;
return [false];
end;;
......@@ -2,13 +2,14 @@ Read("triangulation.gd");
Read("triangulation.gi");
Read("line_intersection.g");
Read("TwoTriangleIntersection.g");
Read("outer_hull.g");
# global variable
eps:=1.*10^-6;
#
calculate_intersections:=function(vof,coordinates,intersections_only)
local l,i,j;
local l,i,j,intersection,intersection2;
l:=[];
for i in [1..Size(vof)] do
l[i]:=[List([1..3],j->coordinates[vof[i][j]]),Set([Set([1,2]),Set([2,3]),Set([3,1])])];
......@@ -97,7 +98,7 @@ end;;
# join the data of all triangles
join_triangles:=function(list)
local vertices,VerticesOfFaces,map,i,pos;
local vertices,VerticesOfFaces,map,i,pos,data;
vertices:=[];
VerticesOfFaces:=[];
for data in list do
......@@ -147,7 +148,7 @@ end;;
# fix_intersections_planar([[[1.,0,0],[4.,0,0],[1.,0,4],[2.,0,0],[1.,0.,3],[2,0,2]],[[1,2],[2,3],[3,1],[4,5],[1,6]]]);
fix_intersections_planar:=function(data)
local i,j,res,entry,entry1,entry2,entry_in_i,entry_in_j,addy,check_edges;
local i,j,res,entry,entry1,entry2,entry_in_i,entry_in_j,addy,check_edges,orig_j,orig_i,cur;
# check for non-paralel intersecting lines
data:=clean_data(data,eps);
check_edges:=Combinations([1..Size(data[2])],2);;
......@@ -158,7 +159,7 @@ fix_intersections_planar:=function(data)
if i <=Size(data[2]) and j<=Size(data[2]) then
res:=LineSegmentIntersection([data[1][data[2][i][1]],data[1][data[2][i][2]]],[data[1][data[2][j][1]],data[1][data[2][j][2]]],eps);
if res[1] then
Print("Edge ",i,": ",data[2][i]," intersects with Edge ",j,": ",data[2][j],"\n");
#Print("Edge ",i,": ",data[2][i]," intersects with Edge ",j,": ",data[2][j],"\n");
#check if res[2] is already in data[1]
entry:=NumericalPosition(data[1],res[2],eps);
if entry <> fail then
......@@ -238,7 +239,7 @@ fix_intersections_planar:=function(data)
if i <=Size(data[2]) and j<=Size(data[2]) then
res:=LineSegmentIntersectionColinear([data[1][data[2][i][1]],data[1][data[2][i][2]]],[data[1][data[2][j][1]],data[1][data[2][j][2]]],eps);
if res[1] then
Print("Edge ",i,": ",data[2][i]," intersects with Edge ",j,": ",data[2][j],"\n");
#Print("Edge ",i,": ",data[2][i]," intersects with Edge ",j,": ",data[2][j],"\n");
#Error();
if res[3]<=4 then
entry1:=NumericalPosition(data[1],res[2][1][2],eps);
......@@ -366,7 +367,7 @@ end;;
#4754
# data contains vertices, edges
my_triangle_fix:=function(data)
local map,i,edge;
local map,i,edge,poss_edges,no_poss_edges,new_edges,ok,e,res;
data:=clean_data(data,eps);
data:=fix_intersections_planar(data);
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment