Intersection de deux sphères

Version imprimable

01/12/2007, 19h25
vinzzzz

Intersection de deux sphères

Hello!

Voila j'ai a ma disposition un ensemble de sphères, et je souhaiterai calculer le volume total qu'elles occupent (je précise que je code en C). Le problème est que ces sphères sont sucéptibles de se chevaucher (sans parler de l'espace vide qui peut apparaitre lorsque deux sphères sont éloignées, mais chaque choses en son temps).

Voila j'ai aucune idée de comment calculer le volume d'intersection de deux sphères (et de ce fait obtenir le volume réellement occupé par ces deux sphères en sommant et en enlevant ce bolume d'intersection).

J'ai fait quelques recherches sur google mais rien pour le moment donc si quelqu'un a une idée (ou un lien) sur la faisabilité de ce truc pour un non spécialiste de la géométrie...

Merci d'avance ;)
01/12/2007, 20h04
Zavonen

Le volume d'une sphère s'obtient en intégrant la surface d'un cercle qui se déplace perpendiculairement à un diamètre.
Ton problème est résolu si tu peux obtenir la surface des deux morceaux de sphère obtenus en coupant la sphère par un plan, puisque l'intersection de deux sphères est exactement la réunion de deux telles 'lentilles'.
Supposons que le plan se trouve à une distance d du centre (inférieure au rayon de la sphère R évidemment).
Tu obtiendras le volume de la lentille en intégrant la surface d'un cercle entre d et R soit la fonction S(x)=pi*(R^2 -x^2).
05/12/2007, 11h24
FR119492

Salut !

Petite remarque concernant la terminologie: une sphère est le lieu des points situés à une distance donnée d'un point appelé centre. C'est donc une surface et l'intersection de deux sphères est un cercle dans l'espace. L'intérieur d'une sphère s'appelle une boule.

Jean-Marc Blanc
05/12/2007, 12h06
pseudocode

Citation:

Envoyé par FR119492

L'intérieur d'une sphère s'appelle une boule.

Et l'interieur d'un cube, ca s'appelle comment ? ;)

Concernant le volume d'intersection sphere/sphere, les formules sont données sur Mathworld:

http://mathworld.wolfram.com/Sphere-...ersection.html
05/12/2007, 12h14
Ulmo

Si c'est juste pour 2 sphères, tu as juste besoin de savoir calculer le volume d'une calotte sphérique (2 calottes en fait). Si ce n'est pas clair j'essayerai de te faire un dessin.
Le volume d'une calotte sphérique est très classique, si tu ne sais pas le faire tu le trouveras par exemple sur Wikipédia.

Par contre si tu regardes l'intersection de 3 sphères ou plus, tu es bien mal parti...

Citation:

Par contre si tu regardes l'intersection de 3 sphères ou plus, tu es bien mal parti...

Dans ce cas une seule solution.
Place toutes tes sphères dans un cube 'maximal'.
Ecris dans un tableau les centres de tes n sphères Cn
Dans un autre tableau les rayons de ces sphères Rn
Le test d'appartenance d'un point à une la sphère n est immédiat (d(M,Cn)<Rn)
Le test d'appartenance à la réunion se fait en itérant le test précédent.
Maintenant recouvre ton cube mettons de côté 1 pour simplifier, avec k^3 points (k à déterminer).
Initialiser V=0;
Faire, par triple boucle, le test d'appartenance à la réunion pour chaque point M(k1,k2,k3), si le point est dans la réunion ajouter 1/k^3 à V.
A la fin V contient une approximation du volume de la réunion, d'autant meilleure que k est grand.
Sinon, pour une grande vitesse avec une précision aléatoire, remplacer cela par Monte-Carlo.

[

Code:

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
#pragma once
 
// le point dans l'espace à 3 dimensions
 
class Point
{public:
double x;
double y;
double z;
public:
	Point(void);
	Point (double,double,double);
	// distance d'un point à un autre
	double distance (Point& );
public:
	virtual ~Point(void);
};
 
#include "StdAfx.h"
#include "Point.h"
#include <math.h>
 
Point::Point(void)
{
}
 
Point::Point(double u,double v,double w)
{x=u;
y=v;
z=w;
}
 
double Point::distance(Point &other)
{return sqrt ((double)((other.x-x)*(other.x-x)+
			  (other.y-y)*(other.y-y)+
			  (other.z-z)*(other.z-z)
			          ));
}
Point::~Point(void)
{
}
 
 
 
#pragma once
#include "point.h"
 
// la classe boule
 
class Boule
{public:
Point centre;
double rayon;
// constructeurs
Boule(void);
Boule(Point&,double);
// coordonnées des points extrêmes
double Topx();
double Botx();
double Topy();
double Boty();
double Topz();
double Botz();
//test d'appartenance
bool contains(Point&);;
~Boule(void);
};
 
#include "StdAfx.h"
#include "Boule.h"
 
Boule::Boule(void)
{
}
Boule::Boule(Point & C, double R)
{centre=C;
rayon=R;
}
bool Boule::contains(Point &P)
{return this->centre.distance(P)<this->rayon;
}
Boule::~Boule(void)
{
}
double Boule::Topx()
{return centre.x+rayon;}
 
double Boule::Botx()
{return centre.x-rayon;}
 
double Boule::Topy()
{return centre.y+rayon;}
 
double Boule::Boty()
{return centre.y-rayon;}
 
double Boule::Topz()
{return centre.z+rayon;}
 
double Boule::Botz()
{return centre.z-rayon;}
 
#pragma once
#include "boule.h"
 
// collection de boules
 
class CollBoules
{
public:
	// le nmbre de boules
	unsigned int nombre;
	// le tableau
	Boule * tabboules;
	// accès à la boule d'indice i (lecture/écriture)
	Boule& operator[](unsigned int );
	//constructeurs
	CollBoules(void);
	CollBoules(unsigned int);
	// destructeur
	virtual ~CollBoules(void);
	// test d'appartenance d'un point
	bool contains(Point&);
	// coordonnées des points extrêmes de la collection
	double Topx();
	double Botx();
	double Topy();
	double Boty();
	double Topz();
	double Botz();
	// calcule le volume de la collection de boules
	// le paramètre k détermine la précision
	// plus k est grand, plus la mesure est précise
	// mais attention ! l'algo est en k^3
	double Volume(unsigned int k);
};
 
#include "StdAfx.h"
#include "CollBoules.h"
 
CollBoules::CollBoules(void)
{
}
CollBoules::CollBoules(unsigned int n)
{nombre=n;
tabboules=new Boule[n];
}
CollBoules::~CollBoules(void)
{delete [] tabboules;
}
Boule& CollBoules::operator[](unsigned int i)
{return this->tabboules[i];
}
bool CollBoules::contains(Point &P)
{for (unsigned int i=0;i<this->nombre;i++)
	if (this->tabboules[i].contains(P))
		return true;
return false;
}
 
double CollBoules::Topx()
{   double result=this->tabboules[0].Topx();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Topx()>result)
			result=this->tabboules[i].Topx();
	return result;
}
 
double CollBoules::Topy()
{   double result=this->tabboules[0].Topy();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Topy()>result)
			result=this->tabboules[i].Topy();
	return result;
}
 
double CollBoules::Topz()
{   double result=this->tabboules[0].Topz();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Topz()>result)
			result=this->tabboules[i].Topz();
	return result;
}
 
double CollBoules::Botx()
{   double result=this->tabboules[0].Botx();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Botx()<result)
			result=this->tabboules[i].Botx();
	return result;
}
 
double CollBoules::Boty()
{   double result=this->tabboules[0].Boty();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Boty()<result)
			result=this->tabboules[i].Boty();
	return result;
}
 
double CollBoules::Botz()
{   double result=this->tabboules[0].Botz();
	for (unsigned int i=0;i<nombre;i++)
		if (this->tabboules[i].Botz()<result)
			result=this->tabboules[i].Botz();
	return result;
}
 
 
double CollBoules::Volume (unsigned int k)
{double V=0;
double deltax=Topx()-Botx();
double deltay=Topy() - Boty();
double deltaz=Topz()-Botz();
double dx=deltax/k;
double dy=deltay/k;
double dz=deltaz/k;
double dv=dx*dy*dz;
Point P;
for (unsigned int i=0;i<k;i++)
for (unsigned int j=0; j<k;j++)
for (unsigned int h=0;h<k;h++)
{P=Point(Botx()+i*dx,Boty()+j*dy,Botz()+h*dz);
if(this->contains(P))V+=dv;
}
return V;
}
 
// programme principal
 
// spheres.cpp : définit le point d'entrée pour l'application console.
//
 
#include "stdafx.h"
#include <stdlib.h>
#include "CollBoules.h"
#include <math.h>
 
int _tmain(int argc, _TCHAR* argv[])
{
	Point P(0,0,0);
	Boule B(P,1);
    CollBoules C(2);
	C[0]=B;
	C[1]=B;
	printf("%lf\n",C.Volume(200));
	system("Pause");
	return 0;
}

06/12/2007, 14h18
Ulmo

Citation:

Envoyé par Zavonen

Le test d'appartenance d'un point à une la sphère n est immédiat (d(M,Cn)<Rn)

C'est un détail, mais il est préférable de comparer les carrés, pour éviter la racine (on stocke le carré des rayons directement, et de toute façon on calcule d'abord le carré de la distance).
06/12/2007, 14h40
Zavonen

Citation:

C'est un détail, mais il est préférable de comparer les carrés, pour éviter la racine (on stocke le carré des rayons directement, et de toute façon on calcule d'abord le carré de la distance).

C'est parfaitement exact ! Il n'y a aucune optimisation.
11/12/2007, 18h46
vinzzzz

OK merci beauocup a tous je vais potasser ca! L'algo donné par Zavonen me semble assez intuitif.

Voici l'implémentation pour Monte-Carlo
Remarque: Pour une boucle à 10000
au lieu de 8000000 pour l'autre méthode
la différence n'est que de 0.03/4

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
double randdouble()
{return ((double) rand())/RAND_MAX;
}
 
double CollBoules::Volume2 (unsigned int k)
{double VP; // volume du parallélépipède
unsigned int N1=0;
double deltax=Topx()-Botx();
double deltay=Topy() - Boty();
double deltaz=Topz()-Botz();
VP=deltax*deltay*deltaz;
double r1,r2,r3;
Point P;
for (unsigned int i=0;i<k;i++)
{  r1=randdouble()*deltax+Botx();
r2=randdouble()*deltay+Boty();
r3=randdouble()*deltaz+Botz();
P=Point(r1,r2,r3);
if(this->contains(P))N1++;
}
return (((double)N1)/k)*VP;
}

12/12/2007, 12h06
vinzzzz

Merci encore! :)