Mageo3d est une bibliothèque de fonctions C++ pour gérer les points et les vecteurs de l'espace euclidien habituel R3. Cette bibliothèque est disponible avec ses sources, l?en-tête : mageo3d.hpp et le développement : mageo3d.cpp, sa version compilée pour Windows : mageo3d.lib, deux programmes de démonstrations : demo1.cpp et demo2.cpp et son manuel d?utilisation : mageo3d.pdf. Le programme demo1.cpp est compilé avec mageo3d.cpp et demo2.cpp avec mageo3d.lib. On a cherché à obtenir un bon style de programmation dans son utilisation. On n?a pas retenu ici l?expression v.scal(w) que l?on trouve souvent sur Internet pour calculer le produit scalaire de 2 vecteurs v et w mais plus simplement : v*w ou scal(v,w). De même, on trouve souvent : v.unit() pour rendre unitaire le vecteur v, mais ici : unit(v) fournit le vecteur unitaire proportionnel à v sans pour autant modifier v. L?avantage principal de Mageo3d est de pouvoir utiliser les fonctions usuelles pour gérer des points et des vecteurs de l?espace euclidien R3 sans avoir à les recréer à chaque fois. On peut aussi assez facilement ajouter, si nécessaire, de nouvelles fonctions dans Mageo3d. Cette nouvelle version de Mageo3d est plus complète et beaucoup plus facile à comprendre, à modifier ou à communiquer.
Source / Exemple :
// ---- fichier mageo3d.hpp -- @author pgl10 ----
#ifndef MAGEO3D_HPP
#define MAGEO3D_HPP
#include <iostream>
#include <cmath>
class V3d {
public:
~V3d();
V3d();
V3d(double a);
V3d(const V3d &v);
V3d(double a, double b, double c);
double getx() const;
double gety() const;
double getz() const;
V3d getxyz() const;
void setx(double a);
void sety(double b);
void setz(double c);
void setxyz(double a, double b, double c);
V3d operator - () const;
V3d operator + (const V3d& v) const;
V3d operator - (const V3d& v) const;
V3d& operator = (const V3d& v);
V3d& operator += (const V3d& v);
V3d& operator -= (const V3d& v);
V3d operator * (double d) const;
double operator * (const V3d &v) const;
V3d& operator *= (double d);
V3d operator / (double d) const;
V3d& operator /= (double d);
V3d operator % (const V3d& w) const;
V3d translate(double dx, double dy, double dz);
V3d translate(const V3d& v);
V3d rotx(double a);
V3d roty(double a);
V3d rotz(double a);
V3d rotate(double a, const V3d& pnt, const V3d& dir);
void print(std::ostream& out) const;
void input(std::istream& in);
private:
double _mx, _my, _mz;
};
V3d unit(const V3d& v);
double norme(const V3d& v);
V3d operator * (double d, const V3d& v);
double scal(V3d& v, V3d& w);
V3d vect(const V3d& v, const V3d& w);
V3d vpts(const V3d& p, const V3d& q);
V3d centre(const V3d& p, const V3d& q);
V3d centre(const V3d& p, const V3d& q, const V3d& r);
V3d centre(const V3d& p, const V3d& q, const V3d& r, const V3d& s);
double dist(const V3d& p, const V3d& q);
double surf(const V3d& p, const V3d& q, const V3d& r);
double volm(const V3d& p, const V3d& q, const V3d& r, const V3d& s);
double angle(const V3d& p, const V3d& q, const V3d& r);
V3d normale(const V3d& p, const V3d& q, const V3d& r);
bool base(V3d b[3], const V3d& v0, const V3d& v1, const V3d& v2);
void duale(V3d a[3], const V3d b[3]);
void comps(double c[3], const V3d& v, const V3d b[3]);
V3d prjd(const V3d& p, const V3d& p1, const V3d& p2);
V3d prjt(const V3d& p, const V3d& p1, const V3d& p2, const V3d& p3);
void pp2d(V3d pq[2], const V3d& p1, const V3d& p2, const V3d& q1, const V3d& q2);
V3d indt(const V3d& p1, const V3d& p2, const V3d& q1, const V3d& q2, const V3d& q3);
void cbary(double c[3], const V3d& p, const V3d& p1, const V3d& p2, const V3d& p3);
std::ostream& operator << (std::ostream& ost, const V3d& v);
std::istream& operator >> (std::istream& ist, V3d& v);
#endif // MAGEO3D_HPP
// ---- fichier mageo3d.cpp -- @author pgl10 ----
// Mageo3d est une bibliothèque de fonctions pour gérer les points ou les vecteurs de R3.
// On obtient un vecteur nul avec les instructions : V3d v; ou : v=V3d(); ou : v=V3d(0);
// On peut aussi faire : v=V3d(1,2,3); ou : w=V3d(v); ainsi que : vx=v.getx()
// et : v.setx(vx); ou même : v.setxyz(a,b,c); On peut aussi faire : t=v; ou : t=-v;
// ou : t=v+w; ou : t=v-w; ou : t*=d; ou : t=v%w pour le produit vectoriel v par w
// L'instruction t=v*w; est valide si v est un double ou un V3d et de même pour w
// t=v/d; est valide aussi. Voir ci-après pour scal(v,w), vect(v,w), p.translate(v)
// et les autres déplacements, voir aussi centre(), dist(), surf(), volm() et normale()
// ainsi que unit(), norme(), angle() et vpts(). Pour créer la base b (V3d b[3]) avec
// les 3 vecteurs v0, v1 et v2 on fera : base(b,v0,v1,v2); et : duale(a,b); pour céer
// sa base duale a, enfin pour calculer les composantes c de v dans b : comps(c,v,b);
// Et pour les points, droites et triangles : prjd(), prjt(), pp2d(), indt() et cbary()
// En mode console, on peut faire : std::cout << v; ou : std::cin >> v;
#include "mageo3d.hpp"
///// Les fonctions membres de la classe V3d /////
V3d::~V3d() {}
V3d::V3d() : _mx(0.0), _my(0.0), _mz(0.0) {}
V3d::V3d(double a) : _mx(a), _my(a), _mz(a) {}
V3d::V3d(const V3d &v) : _mx(v._mx), _my(v._my), _mz(v._mz) {}
V3d::V3d(double a, double b, double c) : _mx(a), _my(b), _mz(c) {}
double V3d::getx() const { return _mx; }
double V3d::gety() const { return _my; }
double V3d::getz() const { return _mz; }
V3d V3d::getxyz() const { return *this; }
void V3d::setx(double a) { _mx = a; }
void V3d::sety(double a) { _my = a; }
void V3d::setz(double a) { _mz = a; }
void V3d::setxyz(double a, double b, double c) { _mx = a; _my = b; _mz = c; }
V3d V3d::operator - () const { return V3d(-_mx, -_my, -_mz); }
V3d V3d::operator + (const V3d& v) const { return V3d(_mx + v._mx, _my + v._my, _mz + v._mz); }
V3d V3d::operator - (const V3d& v) const { return V3d(_mx - v._mx, _my - v._my, _mz - v._mz); }
V3d& V3d::operator = (const V3d& v) { _mx=v._mx; _my=v._my; _mz=v._mz; return *this; }
V3d& V3d::operator += (const V3d& v) { _mx += v._mx; _my += v._my; _mz += v._mz; return *this; }
V3d& V3d::operator -= (const V3d& v) { _mx -= v._mx; _my -= v._my; _mz -= v._mz; return *this; }
V3d V3d::operator * (double d) const { return V3d(_mx*d, _my*d, _mz*d); }
double V3d::operator * (const V3d& v) const { return (_mx*v._mx + _my*v._my + _mz*v._mz); }
V3d& V3d::operator *= (double d) { _mx *= d; _my *= d; _mz *= d; return *this; }
V3d V3d::operator / (double d) const { return V3d(_mx/d, _my/d, _mz/d); }
V3d& V3d::operator /= (double d) { _mx /= d; _my /= d; _mz /= d; return *this; }
V3d V3d::operator % (const V3d& w) const {
double mmx = _my*w._mz - _mz*w._my;
double mmy = _mz*w._mx - _mx*w._mz;
double mmz = _mx*w._my - _my*w._mx;
return V3d(mmx, mmy, mmz);
}
// Attention : Si on déclare : double x, y, z; et : V3d r, v, w;
// l'instruction suivante : r = v + w.translate(x, y, z)
// provoque la translation de w suivie de son addition avec v
// et si on ne veut pas modifier w pendant cette instruction
// il faut faire : r = v + V3d(w).translate(x, y, z);
// Il en est de même pour chacune des transformations suivantes.
V3d V3d::translate(double dx, double dy, double dz) {
_mx = _mx+dx; _my = _my+dy; _mz = _mz+dz;
return *this;
}
V3d V3d::translate(const V3d& v) {
_mx = _mx+v._mx; _my = _my+v._my; _mz = _mz+v._mz;
return *this;
}
V3d V3d::rotx(double a) {
double ca = cos(a);
double sa = sin(a);
double newy = ca*_my - sa*_mz;
double newz = sa*_my + ca*_mz;
_my = newy;
_mz = newz;
return *this;
}
V3d V3d::roty(double a) {
double ca = cos(a);
double sa = sin(a);
double newz = ca*_mx + sa*_mz;
double newx = ca*_mz - sa*_mx;
_mz = newz;
_mx = newx;
return *this;
}
V3d V3d::rotz(double a) {
double ca = cos(a);
double sa = sin(a);
double newx = ca*_mx - sa*_my;
double newy = sa*_mx + ca*_my;
_mx = newx;
_my = newy;
return *this;
}
// pour faire tourner le point concerné d'un angle de t radians
// autour de l'axe défini par le point pnt et le vecteur dir
V3d V3d::rotate(double t, const V3d& pnt, const V3d& dir) {
V3d axe = dir;
double a=axe._mx, b=axe._my, c=axe._mz;
double d=1.0/sqrt(a*a+b*b+c*c);
axe._mx=a*d; axe._my=b*d; axe._mz=c*d;
a=pnt._mx; b=pnt._my; c=pnt._mz;
double u=axe._mx, v=axe._my, w=axe._mz;
double x=_mx, y=_my, z=_mz;
double ct=cos(t), st=sin(t);
_mx=(a*(v*v+w*w)-u*(b*v+c*w-u*x-v*y-w*z))*(1.0-ct)+x*ct+(b*w-c*v-w*y+v*z)*st;
_my=(b*(u*u+w*w)-v*(c*w+a*u-u*x-v*y-w*z))*(1.0-ct)+y*ct+(c*u-a*w-u*z+w*x)*st;
_mz=(c*(u*u+v*v)-w*(a*u+b*v-u*x-v*y-w*z))*(1.0-ct)+z*ct+(a*v-b*u-v*x+u*y)*st;
return *this;
}
void V3d::print(std::ostream& out) const {
out << "[" << _mx << ", " << _my << ", " << _mz << "]";
}
void V3d::input(std::istream& in) {
std::cout << "\nSa composante sur Ox vaut : ";
in >> _mx;
std::cout << "Sa composante sur Oy vaut : ";
in >> _my;
std::cout << "Sa composante sur Oz vaut : ";
in >> _mz;
}
///// Les fonctions simples /////
// unit(v) est le vecteur unitaire proportionnel à v
V3d unit(const V3d& v) {
double vx = v.getx(), vy = v.gety(), vz = v.getz();
double norm2 = vx*vx + vy*vy + vz*vz;
if(norm2 == 0.0) return V3d(0);
double invnrm = 1.0 / sqrt(norm2);
return V3d(vx*invnrm, vy*invnrm, vz*invnrm);
}
// norme(v) est la norme de v
double norme(const V3d& v) {
double vx = v.getx(), vy = v.gety(), vz = v.getz();
double norm2 = vx*vx + vy*vy + vz*vz;
return sqrt(norm2);
}
// d*v est le vecteur qui vaut d fois v
V3d operator * (double d, const V3d& v) { return V3d(d*v.getx(), d*v.gety(), d*v.getz()); }
// produit scalaire de v et w
double scal(V3d& v, V3d& w) {
double vx = v.getx(), vy = v.gety(), vz = v.getz();
double wx = w.getx(), wy = w.gety(), wz = w.getz();
return (vx*wx + vy*wy + vz*wz);
}
// vect(v, w) est le produit vectoriel de v par w
V3d vect(const V3d& v, const V3d& w) {
double vx = v.getx(), vy = v.gety(), vz = v.getz();
double wx = w.getx(), wy = w.gety(), wz = w.getz();
return V3d(vy*wz-vz*wy, vz*wx-vx*wz, vx*wy-vy*wx);
}
// vpts(p, q) est le vecteur (p,q)
V3d vpts(const V3d& p, const V3d& q) {
double px = p.getx(), py = p.gety(), pz = p.getz();
double qx = q.getx(), qy = q.gety(), qz = q.getz();
return V3d(qx-px, qy-py, qz-pz);
}
V3d centre(const V3d& p, const V3d& q) { return 0.5*(p+q); }
V3d centre(const V3d& p, const V3d& q, const V3d& r) { return (p+q+r)/3.0; }
V3d centre(const V3d& p, const V3d& q, const V3d& r, const V3d& s) { return (p+q+r+s)/4.0; }
double dist(const V3d& p, const V3d& q) { return norme(q-p); }
// calcul de la surface du triangle (p,q,r)
double surf(const V3d& p, const V3d& q, const V3d& r) {
double a = norme(r-q), b = norme(p-r), c = norme(q-p);
double d = 0.5*(a+b+c);
return sqrt(d*(d-a)*(d-b)*(d-c));
}
// calcul du volume du tétraèdre (p,q,r,s)
double volm(const V3d& p, const V3d& q, const V3d& r, const V3d& s) {
double base = surf(p, q, r);
V3d norm = normale(p, q, r);
double haut = abs((s-p)*norm);
return base*haut/3.0;
}
// calcul de l'angle en radians entre (p,q) et (p,r)
double angle(const V3d& p, const V3d& q, const V3d& r) {
double a = norme(r-q), b = norme(p-r), c = norme(q-p);
if(b*c==0.0) return 0.0;
double d = (b*b+c*c-a*a)/(2.0*b*c);
return acos(d);
}
// la droite (p1,p2) est-elle valide ?
bool isdr(const V3d& p1, const V3d& p2) {
if(norme(p2-p1)>1.0e-8) return true;
return false;
}
// le triangle (p1,p2,p3) est-il valide ?
bool istr(const V3d& p1, const V3d& p2, const V3d& p3) {
if(norme(p2-p1)<1.0e-8) return false;
if(norme(p3-p2)<1.0e-8) return false;
if(norme(p1-p3)<1.0e-8) return false;
if(sin(angle(p1,p2,p3))<1.0e-8) return false;
return true;
}
// calcul de la normale au triangle (p,q,r)
V3d normale(const V3d& p, const V3d& q, const V3d& r) {
if(!istr(p,q,r)) return V3d(0);
return unit(vect(q-p,r-p));
}
// calcul de la base b avec les 3 vecteurs v0, v1 et v2
bool base(V3d b[3], const V3d& v0, const V3d& v1, const V3d& v2) {
if(abs(vect(v0,v1)*v2)<1.0e-9) {
b[0]=b[1]=b[2]=V3d(0);
return false;
}
b[0]=v0; b[1]=v1; b[2]=v2;
return true;
}
// calcul de la base a duale de la base b
void duale(V3d a[3], const V3d b[3]) {
V3d v;
v = vect(b[1],b[2]); a[0]=v/(b[0]*v);
v = vect(b[2],b[0]); a[1]=v/(b[1]*v);
v = vect(b[0],b[1]); a[2]=v/(b[2]*v);
}
// calcul des composantes c du vecteur v dans la base b
// résultat : v = c[0]*b[0] + c[1]*b[1] + c[2]*b[2]
void comps(double c[3], const V3d& v, const V3d b[3]) {
V3d a[3];
duale(a,b);
c[0] = a[0]*v;
c[1] = a[1]*v;
c[2] = a[2]*v;
}
// calcul de la projection perpendiculaire du point p sur la droite (p1,p2)
V3d prjd(const V3d& p, const V3d& p1, const V3d& p2) {
if(!isdr(p1,p2)) return V3d(0);
V3d v = vpts(p1,p2);
double t = (p-p1)*v/(v*v);
return V3d(p1+t*v);
}
// projection perpendiculaire du point p sur le plan du triangle (p1,p2,p3)
V3d prjt(const V3d& p, const V3d& p1, const V3d& p2, const V3d& p3) {
double c[3];
V3d v0, v1, v2, b[3];
if(!istr(p1,p2,p3)) return V3d(0);
v0 = p2-p1; v1 = p3-p1; v2 = vect(v0,v1);
base(b,v0,v1,v2);
comps(c,p-p1,b);
return V3d(p1+c[0]*v0+c[1]*v1);
}
// calcul de la perpendiculaire commune (p,q) à 2 droites (p1,p2) et (q1,q2)
// résultat : p = pq[0] sur (p1,p2) et q = pq[1] sur (q1,q2)
void pp2d(V3d pq[2], const V3d& p1, const V3d& p2, const V3d& q1, const V3d& q2) {
if(!isdr(p1,p2)) return;
if(!isdr(q1,q2)) return;
V3d vp = p2-p1, vq = q2-q1;
if(norme(vect(vp,vq)) < 1.0e-8*norme(vp)*norme(vq)) return;
// p=p1+dp*p1p2 q=q1+dq*q1q2 pq.p1p2=0 pq*q1q2=0
// d'où : a1+b1*dq+c1*dp=0 a2+b2*dq+c2*dp=0
V3d v1 = q1-p1;
double a1=v1*vp, b1= vq*vp, c1=-vp*vp;
double a2=v1*vq, b2= vq*vq, c2=-vp*vq;
double dp=(a1*b2-b1*a2)/(b1*c2-c1*b2);
double dq=(a1*c2-c1*a2)/(c1*b2-b1*c2);
pq[0]=p1+dp*vp;
pq[1]=q1+dq*vq;
}
// calcul de l'intersection de la droite (p1,p2) avec le plan du triangle (q1,q2,q3)
V3d indt(const V3d& p1, const V3d& p2, const V3d& q1, const V3d& q2, const V3d& q3) {
if(!isdr(p1,p2)) return V3d(0);
if(!istr(q1,q2,q3)) return V3d(0);
V3d pp=p2-p1, nt=normale(q1,q2,q3);
if(abs(pp*nt) < 1.0e-8*norme(pp)) return V3d(0);
double d = ((q1-p1)*nt)/(pp*nt);
return V3d(p1+d*pp);
}
// calcul des composantes barycentriques c de p dans le triangle (p1,p2,p3)
// si p n'est pas dans le plan du triangle (p1,p2,p3) c'est sa projection
// perpendiculaire sur ce plan qui est automatiquement utilisée, résultat :
// p = c[0]*p1 + c[1]*p2 + c[0]*p3 avec c[0]+c[1]+c[2] = 1.0
// quand : c[0]>0 & c[1]>0 & c[2]>0 le point p est dans le triangle (p1,p2,p3)
void cbary(double c[3], const V3d& p, const V3d& p1, const V3d& p2, const V3d& p3) {
if(!istr(p1,p2,p3)) { c[0]=c[1]=c[2]=0.0; return; }
V3d v0 = p2-p1, v1 = p3-p1, v2 = vect(v0,v1);
V3d b[3];
base(b,v0,v1,v2);
double d[3];
comps(d,p-p1,b);
c[0]=1.0-d[0]-d[1];
c[1]=d[0];
c[2]=d[1];
}
// surcharge de : std::cout <<
std::ostream& operator << (std::ostream& ost, const V3d& v) {
v.print(ost);
return ost;
}
// surcharge de : std::cin >>
std::istream& operator >> (std::istream& ist, V3d& v) {
v.input(ist);
return ist;
}
// ---- fichier demo1.cpp -- @author pgl10 ----
#include "mageo3d.hpp"
int main() {
V3d p, q, r, nqr;
std::cout << std::endl << " Demo n\370 1 de Mageo3d" << std::endl << std::endl;
std::cout << " Donnez la position du point P";
std::cin >> p;
std::cout << " Donnez la position du point Q";
std::cin >> q;
std::cout << " Donnez la position du point R";
std::cin >> r;
std::cout << " Le centre du triangle PQR est \205 : " << centre(p,q,r) << std::endl
<< " On peut aussi le retrouver avec : " << (p+q+r)/3.0 << std::endl
<< " La surface du triangle PQR vaut : " << surf(p,q,r) << std::endl
<< " La normale du triangle PQR vaut : " << normale(p,q,r) << std::endl
<< " L'angle entre PQ et PR vaut : " << angle(p,q,r) << " radians" << std::endl;
nqr = unit(vect(normale(p,q,r),vpts(q,r))); // une normale à QR dans le plan de PQR
std::cout << " La hauteur du triangle PQR issue de P vaut : " << abs(scal(vpts(p,q),nqr)) << std::endl;
std::cout << std::endl;
system("pause");
return 0;
}
// ---- fichier demo2.cpp -- @author pgl10 ----
#include "mageo3d.hpp"
#pragma comment(lib, "mageo3d.lib")
int main() {
V3d p, q;
std::cout << std::endl << "Demo n\370 2 de Mageo3d" << std::endl << std::endl;
p = V3d(12, 23, 25);
q = V3d(24, 58, 92);
std::cout << " p : " << p << std::endl
<< " q : " << q << std::endl
<< " p+q : " << p+q << std::endl
<< " q-p : " << q-p << std::endl
<< " p*q : " << p*q << std::endl
<< " p%q : " << p%q << std::endl
<< " 2*q : " << 2*q << std::endl
<< " q*2 : " << q*2 << std::endl
<< " q/2 : " << q/2 << std::endl
<< " vpts(p,q) : " << vpts(p,q) << std::endl
<< " scal(p,q) : " << scal(p,q) << std::endl
<< " vect(p,q) : " << vect(p,q) << std::endl
<< " dist(p,q) : " << dist(p,q) << std::endl
<< " centre(p,q) : " << centre(p,q) << std::endl;
std::cout << std::endl;
system("pause");
return 0;
}
Conclusion :
Toutes les corrections, améliorations ou suppléments sont les bienvenus et toutes les critiques sont intéressantes.
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.