1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
##############################
# Calcule les ou le point d'intersection de deux cercle :
# Cercle 1 de coordonnés (a1,b1) de rayon d1
# Cercle 2 de coordonnés (a2,b2) de rayon d2
# Precision = nombre de chiffre aprés la virgule
##############################
getIntersection=function(a1,b1,d1,a2,b2,d2,precision=5){
# Compute E (cf paper)
computeE=function(a1,b1,a2,b2){
numerateur=-2*a1+2*a2;
denominateur=-2*b1+2*b2;
return(numerateur/denominateur);
}
# Compute F (cf paper)
computeF=function(a1,b1,d1,a2,b2,d2){
numerateur=a1^2+b1^2-a2^2-b2^2-d1^2+d2^2;
denominateur=-2*b1+2*b2;
return(numerateur/denominateur);
}
# Compute intersections if b1 != b2
computeYNotEqual=function(a1,b1,d1,a2,b2,d2){
E=computeE(a1,b1,a2,b2);
F=computeF(a1,b1,d1,a2,b2,d2);
A=E^2+1;
B=(-2*a1+2*E*F+2*E*b1);
C=2*F*b1+b1^2+a1^2+F^2-d1^2;
values=polyroot(c(C,B,A));
x=c();
y=c();
for(i in 1:length(values)){
if(round(Im(values[i]),digit=precision)==0){
x=c(x,Re(values[i]));
y=c(y,-Re(values[i])*E-F);
}
}
if(length(x)==0){
return(NULL);
}
return(matrix(c(x,y),nrow=length(x),ncol=2));
}
# Compute intersections if b1 == b2 and a1 != a2
computeYEqualAndXNotEqual=function(a1,b1,d1,a2,b2,d2){
G=-((a1^2-a2^2-d1^2+d2^2)/(-2*a1+2*a2));
A=1;
B=-2*b1;
C=G^2-2*G*a1+a1^2+b1^2-d1^2;
values=polyroot(c(C,B,A));
x=c();
y=c();
for(i in 1:length(values)){
if(round(Im(values[i]),digit=precision)==0){
x=c(x,G);
y=c(y,Re(values[i]));
}
}
if(length(y)==0){
return(NULL);
}
return(matrix(c(x,y),nrow=length(y),ncol=2));
}
# Compute intersections
if(b1!=b2){
return(computeYNotEqual(a1,b1,d1,a2,b2,d2));
}
else if(a1!=a2){
return(computeYEqualAndXNotEqual(a1,b1,d1,a2,b2,d2));
}
# No intersection found
return(NULL);
}
|
