/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        M A X I M I S A T I O N   D E   L A   P L U S   P E T I T E   D I S T A N C E                                              */
/*        A   L ' I N T E R I E U R   D ' U N   N U A G E   D E   P O I N T S   T R I D I M E N S I O N N E L                        */
/*        E N   C O O R D O N N E E S   C A R T E S I E N N E S                                                                      */
/*        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 "                                              */
/*        A V E C   T E M P E R A T U R E  :                                                                                         */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/recuit_si.43$c' :                                                                                          */
/*                                                                                                                                   */
/*                    Jean-Francois Colonna (LACTAMME, AAAAMMJJhhmmss).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  "INCLUDES.01.I"
                                        /* Introduit le 20051116102340...                                                            */

#define   nPOINTS                                                                                                                       \
                    27                                                                                                                  \
                                        /* Nombre de points du nuage...                                                              */
#define   nITERATIONS                                                                                                                   \
                    90000                                                                                                               \
                                        /* Nombre d'iterations maximum...                                                            */

#define   DELTA_INITIAL                                                                                                                 \
                    80.0
#define   DELTA_FINAL                                                                                                                   \
                    1.0
                                        /* Extrema des variations possibles des coordonnees {x,y}.                                   */
#define   PROBABILITE_DE_PERTURBER_UN_POINT                                                                                             \
                    0.1                                                                                                                 \
                                        /* Probabilite de perturber le point courant...                                              */
#define   PERIODISER_L_UNIVERS                                                                                                          \
                    0                                                                                                                   \
                                        /* L'univers n'est pas periodique et donc limite (0) ou est periodique et donc infini (1).   */
#define   PROBABILITE_D_ACCEPTER_UNE_MAUVAISE_SOLUTION                                                                                  \
                    0.02                                                                                                                \
                                        /* Probabilite d'accepter une mauvaise solution...                                           */

#define   TEMPS                                                                                                                         \
                    (((double)(iterations-1))/((double)(nITERATIONS-1)))                                                                \
                                        /* Simulation d'un temps dans [0,1].                                                         */

extern    double    drand48();
#define   PERTURBATION_ALEATOIRE                                                                                                        \
                    ((drand48() - 0.5)*((delta_initial*(1-temps))+(delta_final*temps)))                                                 \
                                        /* Plus grande variation possible des coordonnees {x,y}, en rappelant que :                  */ \
                                        /*                                                                                           */ \
                                        /*                  drand48() E [0,1]                                                        */ \
                                        /*                                                                                           */ \
                                        /* d'ou la translation de 0.5 qui permet d'avoir un deplacement algebrique dans n'importe    */ \
                                        /* quelle direction et ce de facon isotrope...                                               */

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.                                         */

extern    double    sqrt();
#define   Carre(x)                                                                                                                      \
                    ((x)*(x))                                                                                                           \
                                        /* Definition du carre de 'x'.                                                               */
#define   DistanceEuclidienne(A,B)                                                                                                      \
                    sqrt(Carre(X(B)-X(A)) + Carre(Y(B)-Y(A)) + Carre(Z(B)-Z(A)))                                                        \
                                        /* Definition de la distance entre deux points 'A' et 'B'.                                   */

#define   PERTURBATION(coordonnee,minimum,maximum)                                                                                      \
                    {                                                                                                                   \
                    int       perturber=1;                                                                                              \
                    double    coordonnee_perturbee;                                                                                     \
                    double    temps=TEMPS;                                                                                              \
                                                                                                                                        \
                    while     (perturber == 1)                                                                                          \
                              {                                                                                                         \
                              coordonnee_perturbee = coordonnee + PERTURBATION_ALEATOIRE;                                               \
                                        /* Perturbation de la coordonnee courante...                                                 */ \
                                                                                                                                        \
                              if        (PERIODISER_L_UNIVERS == 0)                                                                     \
                                        {                                                                                               \
                                        }                                                                                               \
                              else                                                                                                      \
                                        {                                                                                               \
                                        if        (coordonnee_perturbee < minimum)                                                      \
                                                  {                                                                                     \
                                                  coordonnee_perturbee = coordonnee_perturbee - minimum + maximum;                      \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                                                                                                                        \
                                        if        (coordonnee_perturbee > maximum)                                                      \
                                                  {                                                                                     \
                                                  coordonnee_perturbee = coordonnee_perturbee - maximum + minimum;                      \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                        }                                                                                               \
                                                                                                                                        \
                              if        ((coordonnee_perturbee < minimum) || (coordonnee_perturbee > maximum))                          \
                                        {                                                                                               \
                                        /* Les coordonnees hors-ecran sont rejetees...                                               */ \
                                        }                                                                                               \
                              else                                                                                                      \
                                        {                                                                                               \
                                        coordonnee = coordonnee_perturbee;                                                              \
                                        perturber=0;                                                                                    \
                                        /* On s'arrete sur la premiere coordonnee perturbee situee dans l'ecran...                   */ \
                                        }                                                                                               \
                              }                                                                                                         \
                    }

#define   DISTANCE_MINIMALE(distance_minimale)                                                                                          \
                    {                                                                                                                   \
                    distance_minimale=+1e+300;                                                                                          \
                                                                                                                                        \
                    for       (i=0 ; i<nPOINTS ; i++)                                                                                   \
                              {                                                                                                         \
                              for       (j=0 ; j<nPOINTS ; j++)                                                                         \
                                        {                                                                                               \
                                        if        (i > j)                                                                               \
                                                  {                                                                                     \
                                                  double    distance_courante=DistanceEuclidienne(Lpoints[i],Lpoints[j]);               \
                                                  distance_minimale=MIN2(distance_minimale,distance_courante);                          \
                                        /* Recherche de la distance minimale courante du nuage de points...                          */ \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                        }                                                                                               \
                              }                                                                                                         \
                    }

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.                                                                 */
static    int       dimZ=0;
#define   Zmin      0
#define   Zmax      (Zmin + (dimZ-1))
                                        /* Definition des profondeurs.                                                               */

#define   IMAGE(x,y)                                                                                                                    \
                    (*(image + (((((int)y)-Ymin)*dimX) + (((int)x)-Xmin))))                                                             \
                                        /* Acces a un point de l'image.                                                              */
#define   PROJECT   0.0
#define   store(n,x,y,z)                                                                                                                \
                    {                                                                                                                   \
                    int       xp=x+((int)(PROJECT*((double)(z))));                                                                      \
                    int       yp=y+((int)(PROJECT*((double)(z))));                                                                      \
                                        /* "Pseudo-projection" du point {x,y,z}.                                                     */ \
                                                                                                                                        \
                    if        (((xp>=Xmin) && (xp<=Xmax)) && ((yp>=Ymin) && (yp<=Ymax)))                                                \
                              {                                                                                                         \
                              IMAGE(xp,yp) = n;                                                                                         \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Rangement d'un point valide d'une image.                                                  */
#define   initialisation(x,y)                                                                                                           \
                    {                                                                                                                   \
                    if        (((x>=Xmin) && (x<=Xmax)) && ((y>=Ymin) && (y<=Ymax)))                                                    \
                              {                                                                                                         \
                              IMAGE(x,y) = NOIR;                                                                                        \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Initialisation de l'image.                                                                */

main()
     {
     int                 generer_les_positions_finales=0;
                                        /* Pour controler la generation des positions finales uniquement (=1) ou des positions       */
                                        /* finales et intermediaires (=0).                                                           */
     int                 sortir_l_image=1;
                                        /* Pour controler la sortie (=1) ou pas (=0) de l'image...                                   */
     int                 x,y;
                                        /* Definition des coordonnees.                                                               */
     unsigned  char      *image;
                                        /* Definition de l'image a generer...                                                        */

     int       i,j,n;
     int       iterations;
                                        /* Index divers...                                                                           */
     point     Lpoints[nPOINTS];
     point     SLpoints[nPOINTS];
                                        /* Liste des points et sa Sauvegarde.                                                        */
     double    delta_initial,delta_final=DELTA_FINAL;
                                        /* Bornes perturbation des coordonnees.                                                      */
     double    distance_minimale_avant;
     double    distance_minimale_apres;
     int       distance_minimale_apres_valide=0;

     Get(dimX,"dimX");
     Get(dimY,"dimY");
     Get(dimZ,"dimZ");
     dimZ = 4*dimZ;
                                        /* Recuperation des dimensions en 'X', en 'Y' et en 'Z' de l'image a generer. On notera la   */
                                        /* multiplication par 4 de 'dimZ' car il est en general 4 fois plus petit que les autres...  */

     image=malloc(dimY*dimX);
                                        /* Definition de l'image a generer...                                                        */

     for  (y=Ymin ; y<=Ymax ; y++)
          {
          for  (x=Xmin ; x<=Xmax ; x++)
               {
               initialisation(x,y);
                                        /* Initialisation de l'image a generer...                                                    */
               }
          }

     delta_initial = DELTA_INITIAL*(((double)MIN2(dimX,dimY))/((double)512));

     for       (n=0 ; n<nPOINTS ; n++)
               {
               double    aveugle1=drand48();
               double    aveugle2=drand48();
               double    aveugle3=drand48();
               double    aveugle4=drand48();
               double    aveugle5=drand48();
               double    aveugle6=drand48();
                                        /* Necessaire (c'est l'experience qui le montre...).                                         */

               X(Lpoints[n]) = Xmin + (0.5*dimX);
               Y(Lpoints[n]) = Ymin + (0.5*dimY);
               Z(Lpoints[n]) = Zmin + (0.5*dimZ);
                                        /* Initialisation arbitraire des points. L'experience montre qu'une configuration            */
                                        /* initiale tres fortement symetrique est preferable a une situation aleatoire...            */
               }

     for       (iterations=1 ; iterations<=nITERATIONS ; iterations++)
               {
               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         X(SLpoints[n]) = X(Lpoints[n]);
                         Y(SLpoints[n]) = Y(Lpoints[n]);
                         Z(SLpoints[n]) = Z(Lpoints[n]);
                                        /* Sauvegarde de l'etat courant des points...                                                */
                         }

               if        (distance_minimale_apres_valide == 0)
                         {
                         if        (iterations == 1)
                                   {
                                   DISTANCE_MINIMALE(distance_minimale_avant);
                                        /* Recherche de la distance minimale avant la perturbation.                                  */
                                   }
                         else
                                   {
                                        /* Dans ce cas 'distance_minimale_avant' est encore valide (calculee lors d'une iteration    */
                                        /* anterieure...).                                                                           */
                                   }
                         }
               else
                         {
                         distance_minimale_avant = distance_minimale_apres;
                         }

               if        (iterations == 1)
                         {
                         fprintf(stderr,"distance minimale initiale=%f\n",distance_minimale_avant);
                         }
               else
                         {
                         }

               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         if        (drand48() < PROBABILITE_DE_PERTURBER_UN_POINT)
                                   {
                                   double    X_perturbe=X(Lpoints[n]);
                                   double    Y_perturbe=Y(Lpoints[n]);
                                   double    Z_perturbe=Z(Lpoints[n]);

                                   PERTURBATION(X_perturbe,Xmin,Xmax);
                                   PERTURBATION(Y_perturbe,Ymin,Ymax);
                                   PERTURBATION(Z_perturbe,Zmin,Zmax);

                                   X(Lpoints[n]) = X_perturbe;
                                   Y(Lpoints[n]) = Y_perturbe;
                                   Z(Lpoints[n]) = Z_perturbe;
                                        /* Perturbation aleatoire des points.                                                        */
                                   }
                         else
                                   {
                                   }
                         }

               DISTANCE_MINIMALE(distance_minimale_apres);
                                        /* Recherche de la distance minimale apres la perturbation.                                  */

               if        (    (distance_minimale_apres < distance_minimale_avant)
                         ||   (drand48() < PROBABILITE_D_ACCEPTER_UNE_MAUVAISE_SOLUTION)
                          )
                                        /* On cherche a maximiser la distance minimale. On notera le test '<' (alors que la          */
                                        /* logique voudrait '<=') a cause du cas ou les points occuperaient initialement tous        */
                                        /* la meme position (d'ou : distance_minimale_avant=0) : dans ce cas, comme un certain       */
                                        /* nombre de points ne sont pas perturbes, on a donc distance_minimale_apres=0, d'ou avec    */
                                        /* avec le test '<=' le rejet de cette solution ; ainsi, on risque d'etre bloque             */
                                        /* indefiniment puisque l'on restaure alors l'etat precedent...                              */
                         {
                         distance_minimale_apres_valide = 0;
                                        /* L'etat anterieur est inchange...                                                          */

                         for       (n=0 ; n<nPOINTS ; n++)
                                   {
                                   X(Lpoints[n]) = X(SLpoints[n]);
                                   Y(Lpoints[n]) = Y(SLpoints[n]);
                                   Z(Lpoints[n]) = Z(SLpoints[n]);
                                        /* Restauration de l'etat anterieur lorsqu'il y a degradation des performances...            */
                                   }
                         }
               else
                         {
                         distance_minimale_apres_valide = 1;
                                        /* On conserve un etat perturbe qui maximise la distance minimale...                         */
                         }

               if        (generer_les_positions_finales == 1)
                         {
                         }
               else
                         {
                         for       (n=0 ; n<nPOINTS ; n++)
                                   {
                                   store(MAX2(NOIR+1,(BLANC*iterations)/nITERATIONS),X(Lpoints[n]),Y(Lpoints[n]),Z(Lpoints[n]));
                                        /* Marquage des points...                                                                    */
                                   }
                         }

               if        (iterations == nITERATIONS)
                         {
                         fprintf(stderr,"distance minimale finale..=%f\n",distance_minimale_avant);
                                        /* ATTENTION : c'est la distance minimale a l'avant-derniere iteration qui est editee car,   */
                                        /* en effet, il n'est pas sur que 'distance_minimale_apres' etait maximale...                */
                         }
               else
                         {
                         }
               }

     if        (generer_les_positions_finales == 1)
               {
               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         store(BLANC,X(Lpoints[n]),Y(Lpoints[n]),Z(Lpoints[n]));
                                        /* Marquage des points...                                                                    */
                         }
               }
     else
               {
               }

     if        (sortir_l_image == 1)
               {
               write(1,image,dimX*dimY);
                                        /* Sortie de l'image...                                                                      */
               }
     else
               {
               }
     }



Copyright © Jean-François Colonna, 2021-2021.
Copyright © CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / Ecole Polytechnique, 2021-2021.