Trouver des bulles de densité dans un nuage de points

**Nebulix** · 21/03/2018, 12h41

A wiwaxia et Sve@r
Je tiens d'abord à redire que j'ai le plus grand respect pour vos personnes et votre travail. Le but de ce forum est cependant d'apporter des réponses originales et d'en débattre avec des avis différents.
Le post original de Sve@r me parait assez représentatif d'un travers trop habituel qui est de ne pas définir assez précisément le but que l'on veut atteindre. Je n'ai toujours pas compris si c'était de mettre une zone en surbrillance sur une image, de donner une liste de disques (x,y,rayon), ou autre.
Je pense que Sve@r aurait gagné beaucoup de temps s'il s'était posé ces questions.
Je maintiens qu'il peut être dangereux de faire des hypothèses sur ce que veut le donneur d'ordre.
Pinaillage ? OUI , je le revendique. En informatique, on ne pinaille jamais assez et çà finit parfois mal.
Faire " un pari sur la compréhension et la curiosité " d'un ordinateur , c'est perdu d'avance.

**wiwaxia** · 22/03/2018, 01h09

Envoyé par Nebulix

... Faire " un pari sur la compréhension et la curiosité " d'un ordinateur , c'est perdu d'avance.

Lorsque l'on donne des leçons de rigueur, il convient pour le moins de s'astreindre à des citations exactes, et non pas caricaturales, afin de ne pas verser dans l'absurdité.
Cette distorsion douteuse dévalue singulièrement ton intervention.

Envoyé par wiwaxia

... Sous ce rapport, toute réponse est un pari sur la compréhension et la curiosité du demandeur, puisqu'il faut passer du langage ordinaire (où les mots ont un sens souvent vague ou multiple) à celui des mathématiques et de la programmation.

Si la question initiale est parfois mal exprimée - ce qui peut naturellement arriver face à un problème difficile, mais qui n'est pas le cas ici - des discussions permettent de mieux cerner les intentions du demandeur, quand elles réunissent des personnes désireuses de s'entendre, et de surmonter les différences langage et de formation.
Dans ces circonstances, il n'est pas interdit (mais au contraire vivement conseillé) de s'y montrer intelligent.

La preuve d'une bonne question, c'est qu'elle reçoit des réponses claires et appropriées; c'est bien ce qui s'est passé jusqu' ici - la seule ombre au tableau étant que l'on attend toujours les tiennes ...

Envoyé par Nebulix

.. En informatique, on ne pinaille jamais assez ...

... et malheureusement ce n'est pas demain la veille qu'on les lira, d'autant que de ton propre aveu - maladresse involontaire

Envoyé par Nebulix

... En informatique, on ne pinaille jamais assez et çà finit parfois mal ...

ta manière de procéder se termine en queue de poisson

.

Il se fait tard; la deuxième solution sera rédigée demain

.

Envoyé par Nebulix

... Je maintiens qu'il peut être dangereux de faire des hypothèses sur ce que veut le donneur d'ordre ...

Il faut savoir vivre dangereusement.

**wiwaxia** · 23/03/2018, 01h11

J'ai repris la distribution des points, initialement concentrée sur des zones trop étroites, et l'échelle des couleurs, insuffisamment contrastée aux fortes concentrations.

On envisage désormais un dénombrement pondéré des (N) points du nuage situés à l'intérieur d'un cercle de rayon donné, et centré en chacun des pixels de l'image.
Il faut pour cela deux nouvelles variables de type tableau, à éléments entiers et dont les dimensions (La, Ha) sont identiques à celles de l'image:
a) un tableau de bytes, initialisé à partir de la liste des points, et constituant une représentation rigoureusement fidèle du nuage; la somme de tous les éléments est égale à (N), et la superposition éventuelle de (k) points au même pixel (x, y) s'y manifeste par une valeur supérieure à l'unité: M[x, y] = k ;
b) un tableau de longints dont les éléments correspondent à la somme coefficientée des points contenus dans le cercle centré en (x, y).
Le plus simple est de rassembler dans une matrice unique les éléments associés par paires, par des déclarations apparentées à ce qui suit:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 CONST DimMax = 2000;
 TYPE PaireBL = RECORD  u: Byte; w: Z_32  END;
      Tab_P   = ARRAY[0..DimMax, 0..DimMax] OF PaireBL;
 VAR Grille: Tab_P;

Le rayon des cercles admettant pour valeur maximale (RdMax), les coefficients de la somme sont par ailleurs calculés et mémorisés dans une matrice carrée de dimensions (2*Rmax + 1); leur valeur ne dépend que de leur distance au centre, ainsi que du rayon du disque, par la relation parabolique:
Coeff[i,j] = R² - i² - j² si (i² + j²) < R²
sinon Coeff[i,j] = 0 ;
d'où les nouvelles déclarations:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 TYPE Tab_C   = ARRAY[-RdMax..RdMax, -RdMax..RdMax] OF Z_32;
 VAR Coeff: Tab_C;

Enfin il y a toujours la liste des points, tableau de vecteurs de dimension (2), dont les composantes sont de type LongInt, et celle des couleurs définie comme tableau de pixels; voici l'ensemble des déclarations:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 CONST NmaxP = 65000; D_Pi = 2 * Pi;
       DimMax = 2000; RdMax = 30; NmaxC = 300;
 
 TYPE Ve2D    = RECORD  x, y: Z_32  END;
      LstVe2D = ARRAY[1..NmaxP] OF Ve2D;
      PaireBL = RECORD  u: Byte; w: Z_32  END;
      Tab_P   = ARRAY[0..DimMax, 0..DimMax] OF PaireBL;
      LstPix  = ARRAY[0..NmaxC] OF Pixel;
      Tab_C   = ARRAY[-RdMax..RdMax, -RdMax..RdMax] OF Z_32;
 
 VAR N_Point, Rdisq, Smax: Z_32; D_Moy: Reel;
     Nuage: LstVe2D;
     Grille: Tab_P;
     ListeC: LstPix;
     Coeff: Tab_C;

La procédure (P2) contient les instructions relatives au calcul de la grille:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 PROCEDURE Aff_G(Np: Z_32);
   CONST C1 = 11; L1 = 19; L2 = L1 + 2; o = 7;
   BEGIN
     E(0015); Wt(C1, L1, 'Nombre de valeurs locales calcul‚es: Nv = ');
     E(0010); Write(Np:7); E(0002); Write((Larg_Image * Haut_Image):(2*o));
 
     E(0015); Wt(C1, L2, 'Valeur maximale:                   Smax = ');
     E(0012); Write(Smax:o); A_
   END;
 
 PROCEDURE Calc_Grille(La, Ha, Rd: Z_32; VAR K_p, K_m: Z_32; VAR G_: Tab_P);
   VAR i, j, Kp, Km, s, x, X1, y, Y1: Z_32; Tx, Ty: Bool;
   BEGIN
     Kp:= 0; Km:= 0;
     FOR x:= 0 TO (La - 1) DO
       FOR y:= 0 TO (Ha - 1) DO
         BEGIN
           s:= 0;               We(69, 19, x, 5);
           FOR i:= -Rd TO Rd DO
             FOR j:= -Rd TO Rd DO
               BEGIN
                 X1:= x + i; Y1:= y + j;
                 Tx:= ((0<=X1) AND (X1<La));
                 Ty:= ((0<=Y1) AND (Y1<Ha));
                 IF (Tx AND Ty) THEN Inc(s, Coeff[i, j] * Grille[X1, Y1].u);
               END;
           Inc(Kp); IF (Km<s) THEN Km:= s; G_[x, y].w:= s
         END;
     K_p:= Kp; K_m:= Km
   END;
 
 PROCEDURE Calc_Coeff(Rd: Z_32; VAR C_: Tab_C);
   VAR i, j, J2, s, t, u: Z_32;
   BEGIN
     FOR i:= -RdMax TO RdMax DO
       BEGIN
         s:= Sqr(Rd); Dec(s, Sqr(i));
         FOR j:= -RdMax TO RdMax DO BEGIN
                                      J2:= Sqr(j); t:= s - J2;
                                      IF (t<0) THEN u:= 0
                                               ELSE u:= t;
                                      C_[i, j]:= u
                                    END
       END
   END;
 
 PROCEDURE Aff_DmRa;
   CONST C1 = 11; C2 = C1 + 42; L1 = 15; L2 = L1 + 2;
   BEGIN
     F(1, L1 - 1, 80, L1 + 7, 2); E(0015);
     Wt(C1, L1, 'Distance moyenne entre les points: Dmoy = ');
     E(0012); Write(D_Moy:7:3);
     E(0015); Wt(C1, L2, 'Rayon du disque (<= ');
     E(0012); Write(RdMax:2);     E(0015); Write(' ):         Rdisc = ');
     E(0008);
     REPEAT
       GotoXY(C2, L2); ClrEol; Read(Rdisq)
     UNTIL (Rdisq<=RdMax);
     E(0010); We(C2, L2, Rdisq, 7)
   END;
 
 PROCEDURE CalcDm(La, Ha, Np: Z_32; VAR D_: Reel);
   VAR Ac: Word; p, q, r: Reel;
   BEGIN
     p:= La * Ha; q:= p / Np; r:= Sqrt(q); D_:= r;
   END;
 
 PROCEDURE Transf_LG(VAR Nu: LstVe2D);
   VAR k: Word; i, j: Z_32;
   BEGIN
     FOR k:= 1 TO NmaxP DO BEGIN
                             i:= Nu[k].x; j:= Nu[k].y;
                             Inc(Grille[i, j].u)
                           END
   END;
 
 PROCEDURE ZeroG(VAR G_: Tab_P);
   CONST Pzero: PaireBL = (u:0; w:0);
   VAR i, j: Word;
   BEGIN
     FOR i:= 0 TO DimMax DO
       BEGIN
         FOR j:= 0 TO DimMax DO G_[i, j]:= Pzero
       END
   END;
 
 PROCEDURE P2;
   VAR Npix: Z_32;
   BEGIN
     ZeroG(Grille); Transf_LG(Nuage);
     CalcDm(Larg_Image, Haut_Image, N_Point, D_Moy);
     Aff_DmRa;      Calc_Coeff(Rdisq, Coeff);
     Calc_Grille(Larg_Image, Haut_Image, Rdisq, Npix, Smax, Grille);
     Aff_G(Npix)
   END;

La procédure (P3) le calcul de la nouvelle image;

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 PROCEDURE Transf_GrIm(La, Ha, Sm: Z_32; VAR Ma: Tab_Pix; VAR G_: Tab_P);
   VAR Xm, Ym: Z_32; k, z: Z_32; f: Reel;
   BEGIN
     f:= NmaxC / Sm;
     FOR Xm:= 0 TO (La - 1) DO
       FOR Ym:= 0 TO (Ha - 1) DO BEGIN
                                   z:= G_[Xm, Ym].w;
                                   k:= Round(f * z);
                                   Ma[Xm, Ym]:= ListeC[k]
                                  END
   END;
 
 PROCEDURE Calc_Palette(VAR L_C: LstPix);
   CONST m = 255;
   VAR Ib, Ir, Iv: Byte; i, j, N1, N2, N3, N4: Word; q: Reel; Px: Pixel;
   BEGIN
     N1:= NmaxC DIV 5; N2:= 2 * N1; N3:= 3 * N1; N4:= 4 * N1;
     FOR i:= 0 TO NmaxC DO
       BEGIN
         q:= i / N1; j:= Trunc(q);
         CASE j OF  0: Ib:= Round(m * q);
                    1: Ib:= m;
                    2: Ib:= Round(m * (3 - q))
                    ELSE Ib:= 0  END;
         CASE j OF  1:    Iv:= Round(m * (q - 1));
                    2..3: Iv:= m;
                    4:    Iv:= Round(m * (5 - q));
                    ELSE  Iv:= 0  END;
         CASE j OF  0..2: Ir:= 0;
                    3:    Ir:= Round(m * (q - 3));
                    ELSE  Ir:= m  END;
         Px[1]:= Ir; Px[2]:= Iv; Px[3]:= Ib; L_C[i]:= Px
       END
   END;
 
  PROCEDURE P3;                                    // Unit‚ Bmp_24_03;
   BEGIN
     F(1, 1, 80, 12, 0); GraphBmp_1(0);
     Calc_Palette(ListeC);
     Transf_GrIm(Larg_Image, Haut_Image, Smax, Matr_Image, Grille);
     Creation_F(Larg_Image, Haut_Image)
   END;

Voici le nuage initial de points (N = 44000),

Nom : 44000_700x700_0.700Nu1x1.png
Affichages : 313
Taille : 49,6 Ko

sa représentation (pour comparaison) par quadrillage d'arête (a = 10), et trois cartes densimétriques correspondant à des disques de rayons respectifs r = 5, 10 et 20):

Nom : 44000_D=3.337_A=10_M=70.png
Affichages : 229
Taille : 20,3 Ko

Nom : 44000_D=3.337_A=10_M=70.png
Affichages : 229
Taille : 20,3 Ko

_

Nom : 44000_D=3.337_R=05__M=000924.png
Affichages : 261
Taille : 98,1 Ko

Nom : 44000_D=3.337_R=10_M=010209.png
Affichages : 220
Taille : 78,9 Ko

_

Nom : 44000_D=3.337_R=20_M=136132.png
Affichages : 246
Taille : 58,3 Ko

Le second procédé supprime la pixelisation de l'image mais pas sa granularité, qui est liée aux fluctuations de la distribution des points.
Il suffit d'insérer deux instructions supplémentaires dans de calcul des couleurs:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 PROCEDURE Transf_GrIm(La, Ha, Sm: Z_32; VAR Ma: Tab_Pix; VAR G_: Tab_P);
   VAR Xm, Ym: Z_32; k, z: Z_32; f: Reel;
   BEGIN
     f:= NmaxC / Sm;
     FOR Xm:= 0 TO (La - 1) DO
       FOR Ym:= 0 TO (Ha - 1) DO BEGIN
                                   z:= G_[Xm, Ym].w;
                                   k:= Round(f * z);
                                   z:= Trunc(0.0199999 * k);
                                   k:= Trunc(60 * z);
                                   Ma[Xm, Ym]:= ListeC[k]
                                  END
   END;

pour faire apparaître des zones de densité bien délimitées (ici r = 20):

Nom : 44000_D=3.337_R=20_M=136132Mult0.0199Mult60.png
Affichages : 246
Taille : 18,5 Ko

Nom : 44000_D=3.337_R=20_M=136132Mult0.0199Mult60.png
Affichages : 246
Taille : 18,5 Ko

Il convient à ce stade de modifier la correspondance effectif/couleur puisque le nombre de teintes employées tombe à 6; comme il est prudent de conserver la palette, on pourra reprendre la dernière procédure de la manière suivante (sous réserve de vérification):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 PROCEDURE Transf_GrIm(La, Ha, Sm: Z_32; VAR Ma: Tab_Pix; VAR G_: Tab_P);
   CONST N_Coul = 6; PlageC = NmaxC DIV (N_Coul - 1);
   VAR Xm, Ym: Z_32; k, z: Z_32; f: Reel;
   BEGIN
     f:= N_Coul / Sm;
     FOR Xm:= 0 TO (La - 1) DO
       FOR Ym:= 0 TO (Ha - 1) DO BEGIN
                                   z:= G_[Xm, Ym].w;
                                   k:= Trunc(f * z);       // k = Trunc(N_Coul * z / Sm) varie de 0 à 6 (pour z = Sm)
                                   IF (k>5) THEN k:= 5;    // Ecrêtage de l'entier (k) à la valeur (N_Coul - 1) = 5
                                   Ma[Xm, Ym]:= ListeC[k * PlageC]
                                  END
   END;

**wiwaxia** · 23/03/2018, 10h11

Le grand espace mémoire occupé par les variables de ce programme m'a d'abord contraint de recourir à la programmation dynamique, en raison de données numériques initialement plus élevées (Structure too large !).
Cela n'avait pas entraîné de grandes complications, la déclaration des variables concernées comportant un tableau de pointeurs:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 CONST DimMax = 2000;
 
 TYPE PaireBL = RECORD  u: Byte; w: Z_32  END;    
      Ligne   = ARRAY[0..DimMax] OF PaireBL;
      P_Ligne = ^Ligne;
      Lst_P   = ARRAY[0..DimMax] OF P_Ligne;
 
 VAR Grille: Lst_P;

L'entier long présent dans la case (x, y) de la grille était donc appelé par Grille[x]^[y].w .

La version la plus simple a été choisie dans le message précédent, pour qu'on s'en tienne à l'essentiel. Chacun fera les transpositions nécessaires dans son propre langage.

Le temps d'exécution de ce nouveau programme est, comme on pouvait s'y attendre, nettement plus important; pour un nuage de 44000 points, le délai de création des images bitmap reste de l'ordre de 5 secondes, tandis que le calcul de la grille apparaît beaucoup plus long, de l'ordre de 1 à 3 minutes pour les valeurs utilisées: 105, 125 et 165 s lorsque le rayon prend les valeurs: 25, 30 et 40 .
Il y a ici un facteur sévèrement contraignant pour qui souhaite un résultat rapide.

De plus, l'intervention de rayons croissants épaissit les contours mais atténue leurs irrégularités; les images ci-dessous correspondent à des disques de rayons 25, 30 et 40 :

Nom : 44000_D=3.337_R=25_M=317152_6Coul.png
Affichages : 213
Taille : 6,9 Ko

Nom : 44000_D=3.337_R=25_M=317152_6Coul.png
Affichages : 213
Taille : 6,9 Ko

_

Nom : 44000_D=3.337_R=30_M=609951_6Coul.png
Affichages : 217
Taille : 7,0 Ko

_

Nom : 44000_D=3.337_R=40_M=1644614_6Coul.png
Affichages : 230
Taille : 7,6 Ko

**Nebulix** · 23/03/2018, 16h50

Nom : c1.JPG
Affichages : 205
Taille : 137,4 Ko

Nom : c2.JPG
Affichages : 237
Taille : 130,0 Ko

Le calcul prend un peu mois de 0.8s, pour des images 2000*2000.(affichées en 500*500).

**Nebulix** · 23/03/2018, 16h52

Erratum : images 1000*1000.
Désolé

**Nebulix** · 23/03/2018, 17h10

Nom : c3.JPG
Affichages : 221
Taille : 128,9 Ko

On peut descendre à 0.2 s en remarquant qu'il ne sert à rien d'avoir une résolution d'image beaucoup plus petite que le "rayon" de la "zone", i.e. en calculant une image finale de 500 * 500.

**wiwaxia** · 23/03/2018, 18h12

Pas mal, l'image obtenue, surtout vis-à-vis du temps de calcul

Combien y a-t-il de points ? Comment le "rayon" de la "zone" est-il défini ? Il est intéressant de connaître les caractéristiques du nuage.

Et surtout quel logiciel as-tu utilisé ? L'information concerne directement l'auteur du forum.

**Sve@r** · 24/03/2018, 10h26

Envoyé par wiwaxia

Et surtout quel logiciel as-tu utilisé ? L'information concerne directement l'auteur du forum.

Je crois qu'il utilise windev. Mais le rendu graphique (ici très joli) n'est pas vraiment important. Moi, par exemple, je serai sous qgis (la librairie, pas le logiciel) et donc j'utiliserai des widgets Qt...

**wiwaxia** · 24/03/2018, 16h28

En raison d'un penchant irrépressible pour les complications inutiles

, je me suis aperçu seulement ce matin qu'une solution beaucoup plus simple et rapide était envisageable: transformer chacun des points du nuage en un halo circulaire, de rayon égal à celui du disque considéré, en incrémentant les valeurs du voisinage d'une quantité égale au coefficient correspondant, dans la matrice d'entiers associée au corps de l'image.

La variable est ici un tableau d'entiers au format LongInt, qui après mise à zéro de tous ses éléments est initialisée par une procédure unique; les déclarations et instructions spécifiques du programme sont les suivantes:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 CONST DimMax = 2000; ... / ...
 
 TYPE Ligne   = ARRAY[0..DimMax] OF Z_32;
      P_Ligne = ^Ligne;
      Lst_P   = ARRAY[0..DimMax] OF P_Ligne;
      Tab_C   = ARRAY[-RdMax..RdMax, -RdMax..RdMax] OF Z_32;  
 ... / ...
 VAR Grille: Lst_P;
     Coeff: Tab_C;
 ... / ...
 PROCEDURE Calc_G(La, Ha, Np, Rd: Z_32; VAR Nu: LstVe2D; VAR Gm: Z_32);
   VAR i, j, k, Xg, Yg: Z_32;
   BEGIN
     FOR k:= 1 TO Np DO
       FOR i:= -Rd TO Rd DO
         FOR j:= -Rd TO Rd DO
           BEGIN
             Xg:= Nu[k].x + i;
             Yg:= Nu[k].y + j;
             IF (((0<=Xg) AND (Xg<La)) AND ((0<=Yg) AND (Yg<Ha))) THEN
               Inc(Grille[Xg]^[Yg], Coeff[i, j])
           END;
     k:= 0;
     FOR Xg:= 0 TO (La - 1) DO
       FOR Yg:= 0 TO (Ha - 1) DO
         IF (k<Grille[Xg]^[Yg]) THEN k:= Grille[Xg]^[Yg];
     Gm:= k
   END;
 
 PROCEDURE ZeroG;
   VAR i, j: Word;
   BEGIN
     FOR i:= 0 TO DimMax DO
       BEGIN
         New(Grille[i]);
         FOR j:= 0 TO DimMax DO Grille[i]^[j]:= 0
       END
   END;
 ... / ...
 PROCEDURE P2;
   VAR Npix: Z_32;
   BEGIN
     ZeroG;    ... / ... Calc_Coeff(Rdisq, Coeff);
     Calc_G(Larg_Image, Haut_Image, N_Point, Rdisq, Nuage, Gmax);
   ... / ...
   END;

Le calcul des coefficients attachés au disque a été légèrement modifié, afin que les termes extrêmes ((Rd, 0), (0, Rd), ... etc) ne soient pas nuls:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 PROCEDURE Calc_Coeff(Rd: Z_32; VAR C_: Tab_C);
   VAR i, j, J2, s, t, u: Z_32;
   BEGIN
     FOR i:= -RdMax TO RdMax DO
       BEGIN
         s:= Round(Sqr(Rd + 0.5)); Dec(s, Sqr(i));
         FOR j:= -RdMax TO RdMax DO BEGIN
                                      J2:= Sqr(j); t:= s - J2;
                                      IF (t<0) THEN u:= 0
                                               ELSE u:= t;
                                      C_[i, j]:= u
                                    END
       END
   END;

Les images ci-dessous représentent un nuage de 44000 points, sur un domaine initial de 700x700 pixels.
Les calculs impliquant des disques de rayons (10, 20, 30) ont été réalisés sur des délais allant de une à cinq secondes: ce n'est pas instantané, mais nettement plus rapide que l'exécution de la version précédente

.

Nom : Sp_44000_0.999_0.150.png
Affichages : 194
Taille : 52,4 Ko

Nom : Sp_44000_D=3.337_R=10_M=013633.png
Affichages : 197
Taille : 8,9 Ko

_

Nom : Sp_44000_D=3.337_R=20_M=112600.png
Affichages : 219
Taille : 11,3 Ko

Nom : Sp_44000_D=3.337_R=30_M=406421.png
Affichages : 177
Taille : 12,6 Ko

_

Nom : Sp_44000_R=20_M=112600.png
Affichages : 211
Taille : 90,5 Ko

Pour la dernière image, en coloration continue, le rayon du disque vaut 20.

Il y a probablement plein de détails à reprendre dans ce programme, compte tenu de toutes les variantes essayées.

**Nebulix** · 26/03/2018, 11h10

Bonjour
Des considérations générales paraissent incongrues quand elles ne sont pas illustrées par des exemples.

Dans ce cas précis, je m'étonne encore que Sver n'ait pas donné d'exemple de nuage à traiter, qu'il ait fallu attendre le 9e post pour préciser :
" En fait je bosse sur un shape géographique...", encore quelques uns pour préciser ce que cela veut dire.

Dans un élan de générosité , wiwaxia et tbc ont "réinventé" la question et proposé des solutions. Il semble qu'ils sont tombés juste, au moins sur le but recherché, de simplement montrer une image destinée à être vue par un observateur.( Pas de donner les barycentres et rayons de giration des zones ...)
En fait, je réponds au problème de wiwaxia, sans vraiment savoir si c'est celui de sver.

Si le but est l'oeil humain, et qu'on veut travailler avec des zones assez grandes, alors il devient évident qu'on n'a pas besoin de calculer beaucoup de pixels, ce qui permet d'accélérer considérablement le calcul.
Si le 'rayon' est 10, le diamètre est 20 et 4 ou 5 pixels suffisent au rendu, ce qui divise le temps de calcul par 16 ou 25 ...
Dans ce cas précis on peut aussi remarquer qu'il vaut mieux faire la boucle principale sur la quantité 'rare', ici les points donnés au départ plutôt que les pixels de l'image. (le dernier post de wiwaxia ?)

J'utilise un logiciel que beaucoup considèrent comme une vieillerie : Delphi 7,( dernière version disponible gratuitement). Le code ( écrit rapidement et pas spécialement optimisé) contient :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
type b3=array[1..3]of byte;
b3a=array   [0..2000]of b3;
pb3a=^b3a;
im=array[0..1000,0..1000]of single;

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
var
bb:b3;
p:pb3a;
bm:tbitmap;
ker:array[-20..20,-20..20]of single; r,r2:single;
ker2:array[-10..10,-10..10]of single;
im2,im1:im;

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
//définition du rayon et de la fonction associée
r2:=sqr(5);
for i:=-10 to 10 do for j:=-10 to 10 do begin
ker2[i,j]:=1/exp(((i*i)+(j*j))/r2);
end;

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
procedure TForm1.Button4Click(Sender: TObject);
var i,j,k,x,y:integer;  zm:single;
begin
  for x:=0 to 1000 do  for y:=0 to 1000 do   im2[x,y]:=0;
  //boucle principale
  for x:=0 to 1000 do  for y:=0 to 1000 do   if im1[x,y]>0 then
  begin
  for i:=-10 to 10 do for j:=-10 to 10 do
  if( ((x div 2+i)>=0) and  ((x div 2+i)<=1000) and ((y div 2+j)>=0) and  ((y div 2+j)<=1000)) then
  im2[x div 2+i,y div  2+j]:=im2[x div 2+i,y div 2+j]+ker2[i,j];
  end;
  //normalisation
   zm:=0;
    for x:=10 to 490 do
  for y:=10 to 490 do
   if zm< im2[x,y]  then  zm:=im2[x,y];
  for x:=10 to 490 do
  for y:=10 to 490 do
    im2[x,y]:=  im2[x,y]/zm*254;
    //copie dans le bitmap
    bm.Height:=500;bm.Width:=500;
       for y := 0 to Bm.Height -1 do
    begin
      P := Bm.ScanLine[y];
      for x := 0 to Bm.Width -1 do begin
      k:=trunc(0.1+im2[x,y]);
//table des couleurs sommaire
if k>200 then begin      bb[1]:=255;  bb[2]:=255; bb[3]:=255; end;//1 bleu  2 vert    3 rouge
if (k<=200)and(k>150)then begin bb[1]:=0;  bb[2]:=255; bb[3]:=0;end;
if (k<=150)and(k>100)then begin bb[1]:=255;  bb[2]:=0; bb[3]:=0;end;
if (k<=100)and(k>50)then begin bb[1]:=0;  bb[2]:=0; bb[3]:=255;end;
if (k<=50)then begin bb[1]:=0;  bb[2]:=0; bb[3]:=0;end;
        P[x] := bb;
    end;   end;
//affichage
with image1 do Canvas.stretchDraw(Rect(0,0, image1.Width-1,image1.Height-1),Bm);
 
end;

Quitte à passer encore pour un pinailleur, je répète qu'être précis dans l'énoncé d'un problème est essentiel à sa solution.

**wiwaxia** · 26/03/2018, 12h51

Envoyé par Nebulix

... wiwaxia et tbc ont "réinventé" la question et proposé des solutions. Il semble qu'ils sont tombés juste, au moins sur le but recherché, de simplement montrer une image destinée à être vue par un observateur.( Pas de donner les barycentres et rayons de giration des zones ...

Il m'a paru évident, dès le premier message, qu'il fallait représenter les variations de la densité moyenne locale définie soit sur une case (représentation en mosaïque), soit sur un disque (représentation polychrome, dégradée ou discontinue); le but étant de repérer la (ou les) région(s) de concentration maximale.

Envoyé par Nebulix

... Dans ce cas précis on peut aussi remarquer qu'il vaut mieux faire la boucle principale sur la quantité 'rare', ici les points donnés au départ plutôt que les pixels de l'image. (le dernier post de wiwaxia ?) ...

Toutes les images sont directement déduites de la liste initiale des coordonnées; seuls varient les paramètres du calcul: arête des cases, ou rayon des disques. Je croyais avoir été suffisamment clair sur ce point.

Envoyé par Nebulix

... J'utilise un logiciel que beaucoup considèrent comme une vieillerie : Delphi 7 ...

C'est encore pire en ce qui me concerne (Pascal), mais c'est parfaitement assumé

.

**wiwaxia** · 27/03/2018, 10h39

@ Nebulix
J'ai lu attentivement les extraits de ton programme, dont la compréhension m'est apparue assez difficile en raison d'une indentation anarchique, et de l'absence totale de constantes (l'une des valeurs numériques y intervient douze fois !).
En utilisant à fond la fonction Search de l'éditeur Virtual Pascal, la mule que je suis a donc repris le texte source, en conservant tout ce qui pouvait l'être, afin d'avoir une idée de l'algorithme sous-jacent ... sans garantie de bon fonctionnement, faute d'avoir pu le compiler.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 CONST Dim_Im = 1000;
       N1 = 10; N2 = 500; N3 = N2 - N1;   // 500 - 10 = 490
       Max = 255;
 var bb:b3;
     p:pb3a;
     bm:tbitmap;
     ker:array[-20..20,-20..20]of single; r,r2:single;     //ker non utilis‚
     ker2:array[-N1..N1,-N1..N1]of single;
     im2,im1:im;
 
//définition du rayon et de la fonction associée
   r2:=sqr(5);
   for i:=-N1 to N1 do
     for j:=-N1 to N1 do begin
                           ker2[i,j]:=1/exp(((i*i)+(j*j))/r2);
                         end;
 
 procedure TForm1.Button4Click(Sender: TObject);
   var i, j, k, K1, x, y:integer; zm:single; TestX, TestY: Boolean;
   begin
     for x:=0 to Dim_Im do
       for y:=0 to Dim_Im do im2[x,y]:=0;        // im2
  //boucle principale
     for x:=0 to Dim_Im do
       for y:=0 to Dim_Im do
         if im1[x,y]>0 then  //im1
           for i:=-N1 to N1 do
             for j:=-N1 to N1 do
               BEGIN
                 TestX:= ((x div 2+i)>=0) and ((x div 2+i)<=Dim_Im);
                 TestY:= ((y div 2+j)>=0) and ((y div 2+j)<=Dim_Im);
                 if (TestX AND TestY) then
                   im2[x div 2+i,y div  2+j]:= im2[x div 2+i,y div 2+j] + ker2[i,j];   // im2
               END;
  //normalisation
     zm:=0;
     for x:=N1 to N3 do
       for y:=N1 to N3 do
         if zm< im2[x,y]  then  zm:=im2[x,y];                  // im2
     for x:=N1 to N3 do
       for y:=N1 to N3 do im2[x,y]:=  im2[x,y]/zm*254;
    //copie dans le bitmap
     bm.Height:=N2; bm.Width:=N2;
     for y := 0 to Bm.Height -1 do
       begin
        P := Bm.ScanLine[y];
        for x := 0 to Bm.Width -1 do
          begin
            k:=trunc(0.1+im2[x,y]);                            // im2
//table des couleurs sommaire
            K1:= (k - 1) DIV 50;
 
            CASE K1 OF  -1..0: begin
                                 bb[1]:=0;  bb[2]:=0; bb[3]:=0
                               end;
                        1: begin
                             bb[1]:=0; bb[2]:=0; bb[3]:=Max
                           end;
                        2: begin
                             bb[1]:=Max; bb[2]:=0; bb[3]:=0
                           end;
                        3: begin
                             bb[1]:=0; bb[2]:=Max; bb[3]:=0
                           end
                        ELSE begin
                               bb[1]:=Max;  bb[2]:=Max; bb[3]:=Max
                             end;//1 bleu  2 vert    3 rouge
 
(*
            if k>200 then begin
                            bb[1]:=Max;  bb[2]:=Max; bb[3]:=Max
                          end;//1 bleu  2 vert    3 rouge
            if (k<=200)and(k>150) then begin
                                         bb[1]:=0; bb[2]:=Max; bb[3]:=0
                                       end;
            if (k<=150)and(k>100) then begin
                                         bb[1]:=Max; bb[2]:=0; bb[3]:=0
                                       end;
            if (k<=100)and(k>50)  then begin
                                         bb[1]:=0; bb[2]:=0; bb[3]:=Max
                                       end;
            if (k<=50) then begin
                              bb[1]:=0;  bb[2]:=0; bb[3]:=0
                            end;                                                     *)
            P[x] := bb;
          end
        end;
//affichage
     with image1 do Canvas.stretchDraw(Rect(0,0, image1.Width-1,
                                       image1.Height-1),Bm)
   end;

J'avoue n'avoir rien compris à la genèse de la variable (im2), dont les éléments sont incrémentés après initialisation à zéro d'une façon qui m'échappe totalement:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

im2[x div 2+i,y div  2+j]:= im2[x div 2+i,y div 2+j] + ker2[i,j];

Cela m'intrigue, et j'aurais bien voulu comprendre ...
La variable jumelle (im1), qui n'intervient qu'une seule fois dans l'expression d'un booléen:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

if im1[x,y]>0 then

est sans doute initialisée en un autre endroit du programme.

Je me permet simplement de contester le calcul de la fonction matricielle:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
//définition du rayon et de la fonction associée
   r2:=sqr(5);
   for i:=-N1 to N1 do
     for j:=-N1 to N1 do begin
                           ker2[i,j]:=1/exp(((i*i)+(j*j))/r2);
                         end;

au sujet duquel on peut exprimer un double reproche:
a) la fonction de Gauss choisie s'annule pratiquement en dehors des 9 éléments centraux:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
1.000   0.368   0.018   0.000 ...     // 1er terme = ker2[0, 0]
0.368   0.135   0.007   0.000 ...
0.018   0.007   0.000   0.000 ...
0.000   0.000   0.000   0.000 ...
0...

ce qui réduit à un ou deux pixels le rayon effectif du disque de diffusion: le nombre de termes réellement modifiés est par conséquent très inférieur à celui de ceux qui font l'objet d'un calcul ((2*R₁ + 1)² = 11² = 121);
ne faudrait-il pas faire intervenir un coefficient (A) conduisant à des valeurs raisonnablement plus élevées, en convenant par exemple que:
Exp(-A*r²) = 1/2 pour r = R/2 ?
b) le disque de diffusion impose à la fonction la symétrie de révolution vis-à-vis de l'élément central (0, 0), alors que le calcul fait intervenir une matrice carrée: il faut donc imposer l'annulation des éléments situés en-dehors du disque; la fonction associée pourrait ainsi être définie comme suit:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//définition du rayon et de la fonction associée
TYPE T_Ker = ARRAY[-N1..N1,-N1..N1]of Single;
VAR Ker: T_Ker;
 
PROCEDURE CalcK(VAR K_: T_Ker);
  CONST R1 = 5; R2 = R1 * R1;
  VAR i, j, k: Integer; A, u: Single;
  BEGIN
    u:= Ln(2); A:= (4 * u)/R2;
    FOR i:= -N1 TO N1 DO
      FOR j:= -N1 TO N1 DO
        BEGIN
          k:= R2; Dec(k, Sqr(i)); Dec(k(Sqr(j));
          IF (k<0) THEN u:= 0
                   ELSE u:= Exp(-A * k);
          K_[i, j]:= u
        END
  END;

Le coefficient vaut alors A = 0.1109035 , et les éléments de la matrice:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
1.000   0.895   0.642   0.369   0.170   0.062   0.000   ...     // 1er terme = ker2[0, 0]
0.895   0.801   0.574   0.330   0.152   0.000   0.000
0.642   0.574   0.412   0.237   0.109   0.000   0.000  
0.369   0.330   0.237   0.136   0.062   0.000   0.000
0.170   0.152   0.109   0.062   0.000   0.000   0.000   
0.062   0.000   0.000   0.000   0.000   0.000   0.000
0.000   ...

**Nebulix** · 27/03/2018, 14h24

Euh ...
Ce bout de programme a été écrit assez vite, dans le but d'optimiser un temps de calcul, pas pour un concours de beauté.
Les variables im1 et im2 sont définies de type im = array [0..1000,0..1000] of single. Où est le problème ?
La fonction gaussienne s'annule d'autant plus vite qu'on oublie le dénominateur, qui se trouve être justement le carré du rayon ...
Mais elle s'annule quand même assez vite pour ne pas avoir à considérer la forme carrée de la matrice, surtout au niveau de précision dont on a besoin ici.( et dire qu'on me reproche de pinailler)
Plus sérieusement, j'ai peut être omis de préciser que les éléments de im1 sont nuls, sauf ceux qui correspondent à un point de la liste. A partir de ce tableau, je peux afficher la première image.
L'idée de l'algorithme est que chacun de ces points va se traduire par une zone gaussienne dans im2. Ce qui s'écrit par la sommation qui te pose problème.

**wiwaxia** · 28/03/2018, 10h21

Envoyé par Nebulix

... j'ai peut être omis de préciser que les éléments de im1 sont nuls, sauf ceux qui correspondent à un point de la liste. A partir de ce tableau, je peux afficher la première image ...

OK, j'ai compris l'idée du programme ... Mais alors les pixels superposés comptent pour un seul.

Envoyé par Nebulix

... Ce bout de programme a été écrit assez vite, dans le but d'optimiser un temps de calcul, pas pour un concours de beauté ...

Ce qui est en cause n'est pas l'esthétique (encore qu'un bonne présentation favorise la lecture et une compréhension rapide), mais la maîtrise et l'intelligibilité de l'algorithme: j'ai pour réflexe d'introduire une constante dès qu'une valeur apparaît pour la deuxième fois, même dans "tout bout de programme écrit rapidement", sachant par expérience que des modifications ultérieures sont probables, et aisément réalisables dans ces conditions; et cette précaution ne ralentit en rien l'écriture du texte, au contraire.
Un étudiant qui aurait farci son programme de 12 données numériques identiques se serait vu retourner sa copie avec un zéro pointé,

avec en prime l'engueulade du correcteur.

Envoyé par Nebulix

... La fonction gaussienne s'annule d'autant plus vite qu'on oublie le dénominateur, qui se trouve être justement le carré du rayon ...
Mais elle s'annule quand même assez vite pour ne pas avoir à considérer la forme carrée de la matrice, surtout au niveau de précision dont on a besoin ici.( et dire qu'on me reproche de pinailler) ...

C'est bien ce que j'avais pressenti, encore que ce fût pour de mauvaises raisons, parce que la fin de l'expression m'avait échappé - d'où les valeurs fausses attribuées à ton code; je m'étais gouré, en effet.
Au passage,

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

ker2[i,j]:= exp-(((i*i)+(j*j))/r2)

conduit au même résultat en évitant la division (oui

, je pinaille moi aussi).

Envoyé par Nebulix

... L'idée de l'algorithme est que chacun de ces points va se traduire par une zone gaussienne dans im2. Ce qui s'écrit par la sommation qui te pose problème.

a) Il y a d'abord un problème de vocabulaire: le rayon du disque est la distance au delà de laquelle les coefficients s'annulent; ceux-ci sont consignés dans une matrice de dimension supérieure au diamètre, ce qui permet de faire du rayon un paramètre ajustable.
Dans ton programme, le paramètre en question est la demi-dimension (N₁ = 10) de la matrice Ker2, et le rayon, qui en dépend (R₁ = N₁/2), représente la distance à laquelle on trouve les éléments égaux à Exp(-1), soit 0.367879 .
b) L'approximation que tu fais est injustifiée, tant du point de vue de la simplification de l'algorithme, que de la rapidité de son exécution; il est difficile de soutenir que l'instruction conditionnelle

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

(i<sup>2</sup> + j<sup>2</sup>) < Rmax<sup>2</sup>

est lourde à coder, et elle dispense de surcroît du calcul de l'exponentielle dans plus d'un cas sur cinq (1 - Pi/4 = 21.5 %) !
c) Le dénombrement et la sommation des termes de la matrice, selon qu'ils appartiennent ou non au disque inscrit dans le carré, conduit aux résultats suivants:
N_int = 317 ; N_ext = 124 ; Somme = 441 = 21²
S_int = 77.1395; S_ext = 0.9478 ; Rapport = 0.01229 ~ 1.23 % .
Par conséquent, sommer tout les termes du tableau introduit une erreur systématique par excès de 1.23 %, incompatible avec la précision des flottants au format Single (2^-22 (?) ~ 2E^-7) et celle de l'échelle des couleurs (0.5/255 ~ 2E-3).
Un calcul entrepris doit être mené correctement jusqu'à son terme, si l'on veut pouvoir s'assurer de son exactitude; le "niveau de précision dont a besoin" a toujours bon dos, c'est l'excuse ███████ floue de █████████████ l'à-peu-près et le prélude à des ennuis ultérieurs, lorsque l'on veut analyser les résultats.

Ce qui compte, du reste, c'est moins l'intervention de la fonction de Gauss que la symétrie de révolution, et la vérification de la condition:
F(r) = 0 pour (r >= R_max); on peut très bien se limiter à des fonctions monotones décroissantes plus simples, comme:
F₁(r) = (1 - r²/R_max²) ou F₂(r) = F₁(r)² ; la continuité est alors assurée au voisinage de (r = R_max), à l'ordre (0) ou (1).

**Nebulix** · 28/03/2018, 22h45

Un étudiant qui aurait farci son programme de 12 données numériques identiques se serait vu retourner sa copie avec un zéro pointé, avec en prime l'engueulade du correcteur.
...

C'est l'excuse fumeuse de la négligence ...

Il me semble que la limite du minimum de courtoisie de ce forum est dépassée. je n'ai donc plus rien à dire.

**wiwaxia** · 29/03/2018, 09h17

Désolé que des réparties aient été perçues comme blessantes; ce n'était pas le but.

Je réagissais à la façon dont étaient gérées les données essentielles que sont les dimensions des tableaux.

Trouver des bulles de densité dans un nuage de points

Algorithmes et structures de données

Discussions similaires

Partager

Partager