/*************************************************************************************************************************************/
/* */
/* C A L C U L D E G E O D E S I Q U E S D ' U N E S U R F A C E */
/* P A R U N E M E T H O D E D U T Y P E " R E C U I T S I M U L E " : */
/* */
/* */
/* Author of '$xtc/geodesiques.01$vv$c' : */
/* */
/* Jean-Francois COLONNA (LACTAMME, 20160406105800). */
/* */
/*************************************************************************************************************************************/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N S T R E S G E N E R A L E S : */
/* */
/*************************************************************************************************************************************/
#include "INCLUDES.01.I"
#define FLOT(x) \
((double)(x)) \
/* Conversion flottante. */
typedef struct point {
double x;
double y;
double z;
} point;
/* Structure de definition d'un point... */
#define X(P) \
(P.x)
#define Y(P) \
(P.y)
#define Z(P) \
(P.z)
/* Definition de l'acces aux coordonnees d'un point. */
#define distance(A,B) \
sqrt(EXP2(X(B)-X(A)) + EXP2(Y(B)-Y(A)) + EXP2(Z(B)-Z(A))) \
/* Definition de la distance entre deux points 'A' et 'B'. */
/*************************************************************************************************************************************/
/* */
/* I N D I C A T E U R S D E C O N T R O L E : */
/* */
/*************************************************************************************************************************************/
int EditerMauvaisesSolutions_=FAUX;
int EditerMeilleuresSolutions=FAUX;
int GenererImagePlan_uv=FAUX;
int Lister_XYZ=VRAI;
#define NombreIterations \
1000
#define DemiNombreEchantillons \
9
double SeuilPerturbation=0.10;
double DiviseurInter_______uv=1.5;
double MultiplicateurInter_uv=1.5;
double AmplitudeAleatoire=0.0025;
int LimiteurWhile=1000;
/* Indicateurs de controle et parametres divers... */
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E L ' E C H A N T I L L O N N A G E : */
/* */
/*************************************************************************************************************************************/
#define nITERATIONS \
NombreIterations
#define NOMBRE_DE_CHIFFRES \
((int)(log10((double)nITERATIONS)+1))
/* Nombre d'iterations maximum et nombre de chiffres necessaires a son edition... */
#define dECHANTILLONS \
1
#define aECHANTILLONS \
((2*DemiNombreEchantillons)+1)
#define nECHANTILLONS(NumeroEchantillons) \
((double)(NumeroEchantillons-dECHANTILLONS))
#define mECHANTILLONS \
(nECHANTILLONS(aECHANTILLONS)/2)
/* Nombre d'echantillons menant de {uD,vD} a {uA,vA} (eux compris...) et autres... */
/* */
/* On notera que 'nECHANTILLONS' doit etre impair a cause de 'REPLIEMENT(...)' qui doit */
/* bien posseder un point "milieu" ('mECHANTILLONS'). */
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D U G E N E R A T E U R A L E A T O I R E : */
/* */
/*************************************************************************************************************************************/
#define ALEATOIRE_1 \
drand48()
#define ALEATOIRE_2 \
(2*(ALEATOIRE_1 - 0.5))
/* Valeur aleatoire dans : */
/* */
/* ALEATOIRE_2 E [-1,+1] */
/* */
/* en rappelant que : */
/* */
/* drand48() E [0,+1] */
/* */
#define REPLIEMENT(x,XM,YM) \
COND(FLOT(x) <= FLOT(XM) \
,FLOT(x) \
,(-((FLOT(YM)/FLOT(XM))*FLOT(x)) + (2*FLOT(YM))) \
) \
/* Repliement en forme de "chapeau pointu"... */ \
/* */ \
/* L'equation lineaire utilisee est (le repliement ayant lieu au point {XM,YM}) : */ \
/* */ \
/* y = (YM/XM).x si x <= XM */ \
/* y = -(YM/XM).x + 2.YM si x > XM */ \
/* */ \
/* (avec "M" pour "Milieu", la definition de 'x' etant dans [0,2.XM]). */
#define ALEATOIRE_3(x,XM,YM,uvm,uvM) \
(AmplitudeAleatoire*REPLIEMENT(x,XM,YM)*ALEATOIRE_2*ABSO(uvM-uvm)) \
/* Valeur aleatoire dont les extrema varient suivant le "chapeau pointu" afin de davantage */ \
/* perturber {u,v} au milieu du segment {{uD,vD},{uA,vA}} qu'a ses extremites... */
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E S S U R F A C E S : */
/* */
/*************************************************************************************************************************************/
#define Fx_Plan_1________(u,v) \
(u)
#define Fy_Plan_1________(u,v) \
(v)
#define Fz_Plan_1________(u,v) \
(Zplan)
double Zplan=0;
/* Equation d'un plan perpendiculaire a l'axe 'OZ'... */
#define Fx_Sphere_1______(u,v) \
(Rayon*sin(u)*cos(v))
#define Fy_Sphere_1______(u,v) \
(Rayon*sin(u)*sin(v))
#define Fz_Sphere_1______(u,v) \
(Rayon*cos(u))
double Rayon=1;
/* Equation d'une sphere... */
#define Fx_Hyperboloide_1(u,v) \
(RayonA*cosh(u)*cos(v))
#define Fy_Hyperboloide_1(u,v) \
(RayonB*cosh(u)*sin(v))
#define Fz_Hyperboloide_1(u,v) \
(RayonC*sinh(u))
double RayonA=1;
double RayonB=1;
double RayonC=1;
/* Equation d'un hyperboloide a une nappe. */
double PonderFx[]={0,0,1};
double PonderFy[]={0,0,1};
double PonderFz[]={0,0,1};
#define PONDERATION(ponderation,fonction) \
COND(ponderation == 0 \
,0 \
,(ponderation*fonction) \
) \
/* Ponderation optimisee d'une fonction (non calculee si la ponderation est nulle). */
#define Fx(u,v) \
(PONDERATION(PonderFx[0],Fx_Hyperboloide_1(u,v)) \
+PONDERATION(PonderFx[1],Fx_Plan_1________(u,v)) \
+PONDERATION(PonderFx[2],Fx_Sphere_1______(u,v)) \
)
#define Fy(u,v) \
(PONDERATION(PonderFy[0],Fy_Hyperboloide_1(u,v)) \
+PONDERATION(PonderFy[1],Fy_Plan_1________(u,v)) \
+PONDERATION(PonderFy[2],Fy_Sphere_1______(u,v)) \
)
#define Fz(u,v) \
(PONDERATION(PonderFz[0],Fz_Hyperboloide_1(u,v)) \
+PONDERATION(PonderFz[1],Fz_Plan_1________(u,v)) \
+PONDERATION(PonderFz[2],Fz_Sphere_1______(u,v)) \
)
/* Equation de la surface reellement utilisee... */
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D U P L A N { u , v } D E P A R A M E T R A G E D E S S U R F A C E S : */
/* */
/*************************************************************************************************************************************/
double uD=0.0,uA=3.14159265358979312;
double vD=0.0,vA=3.14159265358979312;
/* Definition du point de Depart et d'Arrivee dans le plan {u,v}. */
double Liste_u[aECHANTILLONS-dECHANTILLONS+1];
double Liste_v[aECHANTILLONS-dECHANTILLONS+1];
/* Points {u,v} menant de {uD,vD} a {uA,vA} (eux compris...). */
double Liste_u_perturbees[aECHANTILLONS-dECHANTILLONS+1];
double Liste_v_perturbees[aECHANTILLONS-dECHANTILLONS+1];
/* Points {u,v} perturbes menant de {uD,vD} a {uA,vA} (eux compris...). */
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E L ' I M A G E D U P L A N { u , v } : */
/* */
/*************************************************************************************************************************************/
extern void *malloc();
static int dimX=0;
#define Xmin 0
#define Xmax (Xmin + (dimX-1))
/* Definition des abscisses. */
static int dimY=0;
#define Ymin 0
#define Ymax (Ymin + (dimY-1))
/* Definition des ordonnees. */
#define IMAGE(x,y) \
(*(ImagePlan_uv + (((y-Ymin)*dimX) + (x-Xmin)))) \
/* Acces a un point de l'image. */
#define store_xy(n,x,y) \
{ \
if (((x>=Xmin) && (x<=Xmax)) && ((y>=Ymin) && (y<=Ymax))) \
/* Test de validation des coordonnees... */ \
{ \
IMAGE(x,y) = n; \
} \
else \
{ \
} \
}
#define store_uv(n,u,v) \
{ \
int x=(ABSO(u-uD)/ABSO(uA-uD))*Xmax; \
int y=(ABSO(v-vD)/ABSO(vA-vD))*Ymax; \
/* Denormalisation des coordonnees {u,v} en coordonnees {X,Y}... */ \
\
store_xy(n,x,y); \
}
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E L ' E D I T I O N D E S C O O R D O N N E E S : */
/* */
/*************************************************************************************************************************************/
#define Nombre_u \
20
#define Nombre_v \
20
#define Rayon_Surface \
0.01
#define ROUGE_Surface \
BLANC
#define VERTE_Surface \
NOIR
#define BLEUE_Surface \
NOIR
#define Rayon_Geodesique \
0.02
#define ROUGE_Geodesique \
BLANC
#define VERTE_Geodesique \
BLANC
#define BLEUE_Geodesique \
BLANC
#define COMMENTAIRES_______(sequence)
COMMENTAIRES_______(
xtc
$xCc $BsystemeB $xtc/geodesiques.01$vv$c
FilSTmpB FGeOdEsIqUe
Std
$xtc/$aPout | $AW ' { print $1 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$COORD_X
$xtc/$aPout | $AW ' { print $2 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$COORD_Y
$xtc/$aPout | $AW ' { print $3 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$COORD_Z
$xtc/$aPout | $AW ' { print $4 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe.RAYON
$xtc/$aPout | $AW ' { print $5 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$ROUGE
$xtc/$aPout | $AW ' { print $6 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$VERTE
$xtc/$aPout | $AW ' { print $7 } ' | $SE -e 's/^.*=//' > $FGeOdEsIqUe$BLEUE
set NPoInTs=`$CA $FGeOdEsIqUe$COORD_X | $WCl`
$xrv/particule.10$X \
npoints=$NPoInTs \
LISTE_X=$FGeOdEsIqUe$COORD_X \
LISTE_Y=$FGeOdEsIqUe$COORD_Y \
LISTE_Z=$FGeOdEsIqUe$COORD_Z \
isoles=FAUX chainer=VRAI \
LISTE_RAYON=$FGeOdEsIqUe.RAYON \
LISTE_ROUGE=$FGeOdEsIqUe$ROUGE \
LISTE_VERTE=$FGeOdEsIqUe$VERTE \
LISTE_BLEUE=$FGeOdEsIqUe$BLEUE \
N_AU_CARRE=FAUX \
ZOOM=0.8 ROTATION_OX=1 \
AXYZ=1 BXYZ=0 \
Lz=1000 \
chiffres=0 R=$FGeOdEsIqUe.IMAGE \
$formatI
v $FGeOdEsIqUe.IMAGE
FilSTmpE FGeOdEsIqUe
)
/*************************************************************************************************************************************/
/* */
/* P R O E D U R E S D I V E R S E S : */
/* */
/*************************************************************************************************************************************/
#define PRINT_RESULTATS(editer,message) \
{ \
if (editer == VRAI) \
{ \
printf("Iterations %0*d/%0*d : ",NOMBRE_DE_CHIFFRES,NumeroIterations,NOMBRE_DE_CHIFFRES,nITERATIONS); \
printf(message \
,LongueurTotalePerturbee \
,LongueurTotaleInitiale_ \
); \
} \
else \
{ \
} \
}
/* Rangement d'un point valide d'une image. */
/*************************************************************************************************************************************/
/* */
/* P R O G R A M M E P R I N C I P A L : */
/* */
/*************************************************************************************************************************************/
main()
{
int NumeroEchantillons;
int NumeroIterations;
/* Index divers... */
double LongueurTotaleInitiale_=0;
double Inter_uv=sqrt(EXP2(uA-uD)+EXP2(vA-vD))/nECHANTILLONS(aECHANTILLONS);
double MinimumInter_uv;
double MaximumInter_uv;
MinimumInter_uv = Inter_uv/DiviseurInter_______uv;
MaximumInter_uv = Inter_uv*MultiplicateurInter_uv;
unsigned char *ImagePlan_uv;
/* Definition de l'image a generer... */
unsigned char niveau1;
unsigned char niveau2;
/* Definition des niveau utiles... */
if ( (GenererImagePlan_uv == VRAI)
&& ( (EditerMauvaisesSolutions_ == VRAI)
|| (EditerMeilleuresSolutions == VRAI)
|| (Lister_XYZ == VRAI)
)
)
{
fprintf(stderr,"La generation de l'image ne peut avoir lieu.\n");
GenererImagePlan_uv = FAUX;
}
else
{
}
/*************************************************************************************************************************************/
/* */
/* I N I T I A L I S A T I O N D E L ' I M A G E D U P L A N { u , v } : */
/* */
/*************************************************************************************************************************************/
if (GenererImagePlan_uv == VRAI)
{
int x,y;
/* Definition des coordonnees. */
Get(dimX,"dimX");
Get(dimY,"dimY");
/* Recuperation des dimensions en 'X' et en 'Y' de l'image a generer. */
Get(niveau1,"GRIS_2");
Get(niveau2,"GRIS_4");
/* Recuperation des niveaux utiles. */
ImagePlan_uv = malloc(dimY*dimX);
/* Definition de l'image a generer... */
for (y=Ymin ; y<=Ymax ; y++)
{
for (x=Xmin ; x<=Xmax ; x++)
{
store_xy(NOIR,x,y);
/* Initialisation de l'image... */
}
}
}
else
{
}
/*************************************************************************************************************************************/
/* */
/* I N I T I A L I S A T I O N D U P R O C E S S U S : */
/* */
/*************************************************************************************************************************************/
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
double u;
double v;
double u_1;
double v_1;
point PointPrecedent,PointCourant__;
if (NumeroEchantillons == dECHANTILLONS)
{
u = uD;
v = vD;
/* Pour etre sur de ne pas avoir ici d'erreurs d'arrondi... */
}
else
{
if (NumeroEchantillons == aECHANTILLONS)
{
u = uA;
v = vA;
/* Pour etre sur de ne pas avoir ici d'erreurs d'arrondi... */
}
else
{
u = uD + SCAL((uA-uD),nECHANTILLONS(aECHANTILLONS),nECHANTILLONS(NumeroEchantillons));
v = vD + SCAL((vA-vD),nECHANTILLONS(aECHANTILLONS),nECHANTILLONS(NumeroEchantillons));
/* Les listes des echantillons {u,v} sont initialisees par defaut grace a une */
/* interpolation lineaire... */
}
}
Liste_u[NumeroEchantillons] = u;
Liste_v[NumeroEchantillons] = v;
u_1=Liste_u[NumeroEchantillons-1];
v_1=Liste_v[NumeroEchantillons-1];
X(PointPrecedent) = Fx(u_1,v_1);
Y(PointPrecedent) = Fy(u_1,v_1);
Z(PointPrecedent) = Fz(u_1,v_1);
X(PointCourant__) = Fx(u,v);
Y(PointCourant__) = Fy(u,v);
Z(PointCourant__) = Fz(u,v);
LongueurTotaleInitiale_ = LongueurTotaleInitiale_ + distance(PointPrecedent,PointCourant__);
}
/*************************************************************************************************************************************/
/* */
/* G E N E R A T I O N - 1 - D E L ' I M A G E D U P L A N { u , v } : */
/* */
/*************************************************************************************************************************************/
if (GenererImagePlan_uv == VRAI)
{
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
store_uv(niveau2,Liste_u[NumeroEchantillons],Liste_v[NumeroEchantillons]);
/* Generation de l'image avant perturbation... */
}
}
else
{
}
/*************************************************************************************************************************************/
/* */
/* P E R T U R B A T I O N A L E A T O I R E E T I T E R A T I V E D U S E G M E N T { u , v } : */
/* */
/*************************************************************************************************************************************/
for (NumeroIterations=1 ; NumeroIterations <= nITERATIONS ; NumeroIterations++)
{
double LongueurTotalePerturbee=0;
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
double u=Liste_u[NumeroEchantillons];
double v=Liste_v[NumeroEchantillons];
double u_perturbe,du=0;
double v_perturbe,dv=0;
if (ALEATOIRE_1 <= SeuilPerturbation)
/* Choix des points {u,v} a perturber... */
{
int ItererPerturbation=VRAI;
int LimiteurWhileCourant=LimiteurWhile;
while ((ItererPerturbation == VRAI) && (LimiteurWhileCourant > 0))
{
du = ALEATOIRE_3(nECHANTILLONS(NumeroEchantillons)
,mECHANTILLONS,mECHANTILLONS
,uD,uA
);
dv = ALEATOIRE_3(nECHANTILLONS(NumeroEchantillons)
,mECHANTILLONS,mECHANTILLONS
,vD,vA
);
u_perturbe = u+du;
v_perturbe = v+dv;
if (NumeroEchantillons > dECHANTILLONS)
{
double u_perturbe_1=Liste_u_perturbees[NumeroEchantillons-1];
double v_perturbe_1=Liste_v_perturbees[NumeroEchantillons-1];
double CourantInter_uv;
CourantInter_uv = sqrt(EXP2(u_perturbe-u_perturbe_1)
+EXP2(v_perturbe-v_perturbe_1)
);
/* Distance de deux points {u,v} successifs dans le plan {u,v}. */
if ( (CourantInter_uv < MinimumInter_uv)
|| (CourantInter_uv > MaximumInter_uv)
)
{
LimiteurWhileCourant--;
/* Afin de ne pas boucler indefiniment... */
}
else
{
ItererPerturbation = FAUX;
/* Ces tests sont destines a faire que les points {u,v} successifs ne soient ni trop */
/* prets, ni trop loins les uns des autres dans le plan {u,v}... */
}
}
else
{
ItererPerturbation = FAUX;
}
}
if ((u_perturbe < MIN2(uD,uA)) || (u_perturbe > MAX2(uD,uA)))
{
u_perturbe = u;
/* Afin de rester dans [uD,uA]... */
}
else
{
}
if ((v_perturbe < MIN2(vD,vA)) || (v_perturbe > MAX2(vD,vA)))
{
v_perturbe = v;
/* Afin de rester dans [vD,vA]... */
}
else
{
}
}
else
{
u_perturbe = u;
v_perturbe = v;
/* Cas des points {u,v} qui ne sont pas perturbes... */
}
Liste_u_perturbees[NumeroEchantillons] = u_perturbe;
Liste_v_perturbees[NumeroEchantillons] = v_perturbe;
}
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
if (NumeroEchantillons > dECHANTILLONS)
{
double u=Liste_u[NumeroEchantillons],du;
double v=Liste_v[NumeroEchantillons],dv;
double u_perturbe=Liste_u_perturbees[NumeroEchantillons];
double v_perturbe=Liste_v_perturbees[NumeroEchantillons];
double u_perturbe_1=Liste_u_perturbees[NumeroEchantillons-1];
double v_perturbe_1=Liste_v_perturbees[NumeroEchantillons-1];
point PointPrecedentPerturbe,PointCourantPerturbe__;
X(PointPrecedentPerturbe) = Fx(u_perturbe_1,v_perturbe_1);
Y(PointPrecedentPerturbe) = Fy(u_perturbe_1,v_perturbe_1);
Z(PointPrecedentPerturbe) = Fz(u_perturbe_1,v_perturbe_1);
X(PointCourantPerturbe__) = Fx(u_perturbe,v_perturbe);
Y(PointCourantPerturbe__) = Fy(u_perturbe,v_perturbe);
Z(PointCourantPerturbe__) = Fz(u_perturbe,v_perturbe);
LongueurTotalePerturbee = LongueurTotalePerturbee
+ distance(PointPrecedentPerturbe,PointCourantPerturbe__);
}
else
{
}
}
if (LongueurTotalePerturbee < LongueurTotaleInitiale_)
{
PRINT_RESULTATS(EditerMeilleuresSolutions
,"La solution est MEILLEURE : Lperturbee=%.14f < Linitiale=%.14f\n"
);
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
Liste_u[NumeroEchantillons] = Liste_u_perturbees[NumeroEchantillons];
Liste_v[NumeroEchantillons] = Liste_v_perturbees[NumeroEchantillons];
}
LongueurTotaleInitiale_ = LongueurTotalePerturbee;
/* La solution perturbee etant meilleure, elle remplace la solution initiale... */
}
else
{
PRINT_RESULTATS(EditerMauvaisesSolutions_
,"La solution est mauvaise. : Lperturbee=%.14f >= Linitiale=%.14f\n"
);
}
}
/*************************************************************************************************************************************/
/* */
/* G E N E R A T I O N - 2 - D E L ' I M A G E D U P L A N { u , v } : */
/* */
/*************************************************************************************************************************************/
if (GenererImagePlan_uv == VRAI)
{
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
store_uv(niveau1,Liste_u[NumeroEchantillons],Liste_v[NumeroEchantillons]);
/* Generation de l'image apres perturbation... */
}
write(1,ImagePlan_uv,dimX*dimY);
/* Sortie de l'image presentant le plan {u,v} avant et apres perturbation... */
}
else
{
}
/*************************************************************************************************************************************/
/* */
/* G E N E R A T I O N - 2 - D E L ' I M A G E D U P L A N { u , v } : */
/* */
/*************************************************************************************************************************************/
if (Lister_XYZ == VRAI)
{
double u;
double v;
for (u=uD ; u<=uA ; u=u+((MAX2(uD,uA)-MIN2(uD,uA))/((double)Nombre_u)))
{
for (v=vD ; v<=vA ; v=v+((MAX2(vD,vA)-MIN2(vD,vA))/((double)Nombre_v)))
{
printf("X=%.14f Y=%.14f Z=%.14f RAYON=%f ROUGE=%d VERTE=%d BLEUE=%d \n"
,Fx(u,v)
,Fy(u,v)
,Fz(u,v)
,Rayon_Surface
,ROUGE_Surface
,VERTE_Surface
,BLEUE_Surface
);
/* Edition des coordonnees {X,Y,Z} de quelques points de la geodesique. */
}
}
for (NumeroEchantillons=dECHANTILLONS ; NumeroEchantillons <= aECHANTILLONS ; NumeroEchantillons++)
{
printf("X=%.14f Y=%.14f Z=%.14f RAYON=%f ROUGE=%d VERTE=%d BLEUE=%d \n"
,Fx(Liste_u[NumeroEchantillons],Liste_v[NumeroEchantillons])
,Fy(Liste_u[NumeroEchantillons],Liste_v[NumeroEchantillons])
,Fz(Liste_u[NumeroEchantillons],Liste_v[NumeroEchantillons])
,Rayon_Geodesique
,ROUGE_Geodesique
,VERTE_Geodesique
,BLEUE_Geodesique
);
/* Edition des coordonnees {X,Y,Z} de quelques points de la geodesique. */
}
}
else
{
}
}