Discussion :

[BBouille]

Bonjour à tous,
Après avoir essayé plusieurs méthodes, régression linéaire, moments,... je ne parviens toujours pas à trouver l'angle moyen d'une image inclinée telle que celle jointe.
Je me trompe sûrement mais on dirait qu'il y a moyen de déterminer la pente des traits de cette image.
Si quelqu'un avait une idée je lui en serais vraiment reconnaissant après le temps y consacré.
Cordialement.
BB

[wiwaxia]

Bonjour,

Tu peux procéder de la façon suivante:
a) faire une image agrandie de la grille, et ouvrir le fichier avec Paint;
b) dessiner le plus long rectangle possible (ABCD), en ajustant le coin supérieur gauche (A) sur le noeud le plus à gauche;
c) tracer la demi-droite issue de (A), se superposant à la ligne correspondante de la grille ; elle coupe le côté droit du rectangle en (K);
d) l'inclinaison (t) est alors donnée par la relation: Tan(t) = BK/AB .

On trouve dans le cas présent:
A = ( 31 , 488)
B = (540 , 488)
K = (540 , 526)
Tan(t) = BK/AB = (526 - 488)/(540 - 31) = 38/509 = 0.0746562
t = 4.27°
Le résultat est d'autant plus précis que l'image est plus large, et que les tracés sont plus fins.

[BBouille]

Bonjour,
Merci mais je n'ai peut-être pas été assez complet, je recherche un algorithme qui permettrait d'isoler les lignes et de trouver leur angle et pas leur faire manuellement.
BB

[wiwaxia]

... je recherche un algorithme qui permettrait d'isoler les lignes et de trouver leur angle ...

Alors il faut indiquer quelle est la source des données: un fichier image ? un tableau des coordonnées des noeuds du réseau ?

[BBouille]

Je pars d'une image telle que celle-là.
Je la transforme en noir et blanc et chaque pixel est aussi sauvegardé dans une matrice (x,y) avec pour valeur 0 pour le noir et 1 pour le blanc. Les lignes peuvent être plus ou moins épaisses.
J'ai utilisé la transformée de Radon et de Hough mais sauf mauvaise interprétation de leur algorithme c'est très lourd en termes de calcul et donc de temps.
C'est pourquoi je cherche une autre méthode pour notamment repérer les intersections, extrémités,...
BB

[wiwaxia]

Tu disposes donc d'une matrice (carrée ?) de dimensions (m×m), dont les éléments sont des booléens ou des bytes, et sont caractérisés par deux indices (x, y) variant de zéro à (m - 1).

1°) Soit (N) le nombre approximatif de carreaux correspondant à une dimension de l'image. On peut envisager de localiser le noeud le plus proche de l'élément central (x_I = y_I = (m - 1) DIV 2), et qui doit se trouver dans le domaine [x_I - L , x_I + L]×[y_I - L , y_I + L]
(en convenant de poser L = Round(0.8*m/N)).
Il suffit de compter les éléments non-nuls, dont la distance à l'élément considéré ne dépasse pas un certain seuil (L') nettement inférieur à (L) compte tenu de l'épaisseur finie mais relativement faible des traits de la grille. Cela revient à calculer les (2L + 1)² éléments entiers d'une seconde matrice (M1), par l'exécution d'un code apparenté à ce qui 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
FOR j:= -L TO L DO
  BEGIN
    x:= XI + j;
    FOR i:= -L TO L DO
      BEGIN
        y:= YI + i; k:= 0;
        FOR Dx:= -L1 TO L1 
          FOR Dy:= -L1 TO L1 DO
            IF M[x + Dx, y + Dy]>0 THEN
              BEGIN
                D2:= Sqr(Dx) + Sqr(Dy); 
                IF (D2<Limite) THEN Inc(k)
              END;
        M1[j, i]:= k
      END  
  END;

L'élément (j₀, i₀) de (M1) de valeur maximale correspond au noeud recherché (x₀, y₀).

2°) Envisager un faisceau de droites, passant par le point précédent, et dont l'inclinaison par rapport à l'axes horizontal varie de -45 à +45°; compter les points les plus proches, situés à l'intérieur d'une bande à bords parallèles, admettant pour médiane la droite considérée (largeur à déterminer); on doit observer un maximum prononcé lorsque la droite est parallèle à l'un des axes de la grille, nettement plus accentué que dans le cas précédent.

[BBouille]

Merci Wiw,
Voilà la méthode qui devrait nous amener au but recherché.
Je comprends le principe mais j'ai quand même deux questions à te poser.
Tout d'abord la matrice peut-être rectangulaire mais ça ce n'est pas grave, j'adapterai ton code à une matrice (m,n).

Tu dis: "Soit (N) le nombre approximatif de carreaux correspondant à une dimension de l'image."
Mais N = m dans ton cas non?
Tu dis: "en convenant de poser L = Round(0.8*m/N)".
Pourquoi choisis-tu 0.8?

Merci encore.
BB

[wiwaxia]

... Tout d'abord la matrice peut-être rectangulaire mais ça ce n'est pas grave, j'adapterai ton code à une matrice (m,n) ...

Je m'en doutais, mais me suis limité à la réponse la plus simple.

... Tu dis: "Soit (N) le nombre approximatif de carreaux correspondant à une dimension de l'image."
Mais N = m dans ton cas non? ...

L'image fournie a pour dimensions 158px × 143px, soit en moyenne: m = 150 pixels; et l'on dénombre sur sa largeur et sa hauteur 9.5 et 9 carreaux, d'où : N ~ 9 .
On dispose ainsi d'une estimation du paramètre de la maille carrée du réseau ponctuel ainsi constitué:

a ~158/9.5 ~ 143/9 ~ 16.3 pix .

... Tu dis: "en convenant de poser L = Round(0.8*m/N)".
Pourquoi choisis-tu 0.8?

Chaque noeud étant associé à un carreau de la grille d'arête (a ~ m/N), on doit statistiquement en trouver un dans de domaine central précédemment défini [x_I - L , x_I + L]×[y_I - L , y_I + L] , si l'on prend L = 0.5*a .
afin que les noeuds détectés n'aient pas le mauvais goût de se situer juste sur la limite du domaine, il faut effectuer la recherche sur un domaine plus large, de dimensions

2*L ~1.6*a , d'où: L ~ 0.8*m/N .

L'aire étant de l'ordre de 2.6, on devrait trouver deux zones de densité maximale correspondant à deux noeuds, dont on ne retiendra que le plus proche du centre.

[BBouille]

Ahh d'accord ben je n'avais pas compris alors
Tu recherches en fait les carreaux dans la grille carrée.
Donc ma grille de départ quadrillée induit en erreur d'autant plus que je n'ai pas été assez clair.
Au départ je pensais qu'on pouvait trouver la pente des droites à partir de projection ou en calculant des moyennes ou que sais-je, un peu comme une régression linéaire.
J'avais déjà fait une méthode proche de la tienne mais sur l'épaisseur des traits.
Je vais donc appliquer ton idée et encore réfléchir sur la méthode mathématique si toutefois elle est possible.
Je reviendrai vers toi pour j'espère te donner un résultat valable dans tous les cas de figure.
Encore merci pour ton aide.
BB

[wiwaxia]

Il y a une variante plus simple, en ce qui concerne la première étape: on peut envisager de déterminer le point d'intersection de l'une des droites de la grille avec une verticale située près du bord gauche de l'image, d'abscisse (x₀); on retiendra le point (A) le plus proche de la médiane horizontale, d'ordonnée (Ha/2).

On envisagera par la suite le faisceau de demi-droites (Au) dont l'inclinaison (t) par rapport à l'horizontale varie de (-45°) à (+45°), et le dénombrement des points situés sur le segment correspondant (AM) ou très près de celui-ci (critère non précisé); un pic très étroit doit caractériser la coïncidence de (Au) avec l'une des droites de la grille (segment AM₀).

Compte tenu de la répartition des couleurs

on peut en pratique recourir à un procédé statistique, en accord avec l'inévitable flou du dessin: pour chacune des lignes d'ordonnée (Ym), on calcule la somme (s) des compléments des composantes (r, v, b) des pixels présents dans les (Nc) premières colonnes de l'image; (Nc) étant un nombre impair, le résultat obtenu correspond (à un facteur près) à la valeur moyenne du pixel d'abscisse x₀ - (Nc - 1)/2 .
L'essentiel du code est donné ci-dessous: il détermine la matrice de pixels (Ma2) du corps de la seconde image à partir de celle de la première (Ma1):

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
CONST Nc = 15; Max = 255; NcMax = Nc * Max; X0 = (Nc - 1) DIV 2;
               P100: Pixel = (Max, 0,0); // Pixel rouge
               P010: Pixel = (0, Max, 0); // Pixel vert
... / ...
M1:= Round(0.45 * (Ha - 1)); M2:= Round(0.55 * (Ha - 1)); H1:= Round(0.95 * (Ha - 1));
... / ...
FOR Ym:= 0 TO (Ha - 1) DO
  BEGIN
    s:= 0;
    FOR Xm:= 0 TO (Nc - 1) DO
      FOR k:= 1 TO 3 DO Inc(s, Max - Ma1[Xm, Ym][k]);  // on incrémente (s) de (255 - r), (255 - v), (255 - b): compléments des composantes
                                                       // du pixel Ma1[Xm, Ym] de la 1<sup>re</sup> image, de coordonnées (Xm, Ym)
    h:= Round(Kh * s);
    FOR y:= 0 TO h DO Ma2[Ym, y]:= P010;
    IF ((M1<=Ym) AND (Ym<=M2)) AND (h>H1) THEN
      BEGIN
        Y0:= Ym;
        FOR y:= 0 TO h DO Ma2[Ym, y]:= P100;
      END
  END;                    // Code à vérifier

Voici ce que l'on obtient pour les premières valeurs impaires du nombre (Nc) de colonnes:

Les 10 pics présents se détachent d'autant plus finement sur le bruit de fond que le nombre de colonnes est plus élevé - jusqu'à 15 .
Celui qui nous intéresse (trait rouge) a pour ordonnée Y0 = 73; on retiendra pour l'abscisse (X0) la valeur moyenne: 6 .

[AbsoluteLogic]

Bonjour,

Je tenais à préciser, si j'ai bien compris la solution de Wiwaxia, que l'on peut l'optimiser en ne considérant que les points M qui se situent sur des pixels noirs. Je dirais même vu l'épaisseur des lignes qu'on considère uniquement le pixel noir du milieu lorsqu'il y a un "groupe" de pixels noirs. Et en ne recherchant ces points M qu'à partir de l'ordonnée égale à celle de A et en s'éloignant petit à petit, on a une chance de ne parcourir qu'1 ou 2 segments [AM] pour trouver notre inclinaison.

[wiwaxia]

Bonjour,

... Je tenais à préciser, si j'ai bien compris la solution de Wiwaxia, que l'on peut l'optimiser en ne considérant que les points M qui se situent sur des pixels noirs ...

On pourrait effectivement reprendre le procédé précédent au niveau de la médiane verticale ou près de celle-ci, pour localiser les maximums (points M₁, M₂, M₃ ...); encore faut-il s'assurer pour lequel de ces (N) points le segment (AM_k) est colinéaire à l'un des axes de la grille ...

... Je dirais même vu l'épaisseur des lignes qu'on considère uniquement le pixel noir du milieu lorsqu'il y a un "groupe" de pixels noirs ...

Cela suppose que le groupe de pixels noirs détectés comporte un nombre impair d'éléments ...
Plus généralement, cette distinction dichotomique (noir/blanc) me paraît inadaptée, dans la mesure où l'image source comporte - outre des imperfections locales - des pixels gris et des lignes crénelées dès que celles-ci sont inclinées:

d'où la nécessité de sommer les indices de couleur sur plusieurs colonnes, ce qui revient à calculer un niveau de gris moyen.
Les graphes précédents ont été donnés à titre exploratoire, et il semble que l'on obtienne des résultats fiables (et des pics à peu près symétriques) pour un nombre de colonnes égal à 5, 7, 9 ou 11.
Il faudrait avoir une vue d"ensemble sur les grilles rencontrées pour exprimer des choix plus précis.

... Et en ne recherchant ces points M qu'à partir de l'ordonnée égale à celle de A et en s'éloignant petit à petit, on a une chance de ne parcourir qu'1 ou 2 segments [AM] pour trouver notre inclinaison.

Sauf si les lignes de la grille présentent des orientations voisines de celles des bissectrices du repère, ce que n'exclut pas l'énoncé du sujet ...

Il faudra alors balayer le secteur [-45° ; +45°] (et même un peu plus !) pour trouver l'inclinaison ...

[Nebulix]

Peux tu faire une transformée de Fourier à 2 D ?

[wiwaxia]

Peux tu faire une transformée de Fourier à 2 D ?

Je crains que cela conduise à l'étude d'une fonction de deux variables F(u, v) donnée par une intégrale double, donc résultant d'une double sommation sur le domaine de l'image ... On risque fort de l'enliser dans les calculs très lourds, ceux-là mêmes que l'auteur du sujet tient à éviter.
Le procédé en deux étapes proposé (recherche d'un point de la grille, puis de la pente de la droite de densité maximale) - laborieux, j'en conviens - vise justement à déterminer les caractéristiques d'une fonction bipériodique.

Peut-être as-tu en tête une solution facile à mettre en oeuvre? Il serait intéressant de la connaître

[BBouille]

Holà Wiw je venais te dire où j'en étais et je vois que tu as encore consacré du temps à mon problème.
Je suis passé par ta méthode en projetant les pixels noirs sur des droites inclinées.
J'avais placé une droite sur le coin supérieur gauche que je faisais pivoter de 0 à 45° et la même chose sur le coin inférieur gauche.
Mais l'inconvénient c'est de nouveau le temps de calcul, il faut balayer plusieurs fois la grille.
Finalement enfin j'espère, j'ai programmé un truc relativement simple, j'ai encore quelques petits soucis de pentes positives ou négatives.
On balaie l'image de haut en bas et de gauche à droite et dès qu'on tombe sur un pixel noir on le mémorise dans une variable de type tableau de Pt(x,y), on recherche le pixel noir adjacent (NE-E-SE) ou (NO-O-SO) qui offre la pente la plus proche de celle de la droite qu'on forme petit à petit.
On met chaque pixel mémorisé à blanc. On propage la ligne ainsi créée jusqu'à ne plus rencontrer de pixel noir.
On repart de l'endroit où on avait rencontré le premier pixel noir et on recommence le processus jusqu'au dernier pixel noir.
À la fin on se retrouve avec un ensemble de droites voire de courbes.
On élimine les droites les plus longues au cas où elles seraient issues d'un bord apparaissant parfois lors du scannage.
On compare la pente de quelques-unes des droites et alors on peut fixer la pente générale de la grille.
De cette façon on peut également reconstituer la grille de façon numérisée.
J'ai toujours l'impression mais sûrement fausse qu'on peut trouver par calcul la pente générale d'une image un peu comme une régression linéaire mais je vais m'arrêter là comme je ne vois pas comment faire mieux qu'en un seul balayage.
Si je ne te lis plus encore merci pour ton dévouement.
BB

[wiwaxia]

... À la fin on se retrouve avec un ensemble de droites voire de courbes.
On élimine les droites les plus longues au cas où elles seraient issues d'un bord apparaissant parfois lors du scannage.
On compare la pente de quelques-unes des droites et alors on peut fixer la pente générale de la grille.
De cette façon on peut également reconstituer la grille de façon numérisée ...

Les idées sont pertinentes, mais leur mise en oeuvre apparaît très incertaine, compte tenu de la diversité des cas rencontrés, et de la nécessité de coder tout traitement manuel.

J'ai toujours l'impression mais sûrement fausse qu'on peut trouver par calcul la pente générale d'une image un peu comme une régression linéaire mais je vais m'arrêter là comme je ne vois pas comment faire mieux qu'en un seul balayage ... BB

La régression linéaire peut conduire à de résultats précis, mais à condition d'avoir préalablement identifié les noeuds théoriquement alignés, donc les droites du réseau ponctuel de densité maximale.
C'est possible à condition de balayer la plus grande partie de l'image (en laissant tomber les 5 ou 10 % des points situés au plus près des bords), et de compter les pixels noirs les plus proches en utilisant une fonction de pondération dépendant de la distance, et s'annulant au-delà d'un certain seuil (d_max) - comme par exemple:

Coef(d) = (1 - (d/d_max)²)⁴ .

La fonction de deux variables F(x, y) représentable sur une image de mêmes dimensions (La, Ha) que celle de la grille pourrait être donnée par un algorithme de la forme suivante (rédigé à la volée , donc à vérifier):

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
VAR D2, Ic, u, v, Xmax, Xmin, Ymax, Ymin: LongInt;
    C, s: Exrended;
Xmin:= Round(0.05*La); Xmax:= Round(0.95*La);
Ymin:= Round(0.05*Ha); Ymax:= Round(0.95*Ha);
 
FOR x:= Xmin TO Xmax DO
  FOR y:= Ymin TO Ymax DO
    BEGIN
      s:= 0;
      FOR u:= -h TO h DO
        FOR v:= -h TO h DO
          BEGIN                               // Ma1: matrice de pixels du corps de la première image
            Ic:= Ma1[x + u, y + v][1];  Inc(Ic, Ma1[x + u, y + v][2]); Inc(Ic, Ma1[x + u, y + v][3]); // Ic = somme des trois indices de couleur
            D2:= Sqr(u) + Sqr(v);        C:=  Coef(D2);  // Calcul du carré de la distance, et du coefficient de pondération
            s:= s + (C * Ic)     
          END;
      Ma2[x, y]:= FonctionCouleur(s)          // Ma2: matrice de pixels du corps de la seconde image
    END;
 
// il faut h <= Xmin et Ymin

La génération de l'image est intéressante, mais il faudra déterminer ensuite l'orientation des axes de la grille, sur lesquels la distance entre deux noeuds consécutifs et minimale ... Tu n'es peut-être pas au bout de tes peines

La démarche la plus simple me paraît être la suivante:
a) calcul de l'ordonnée moyenne d'une droite pour une valeur donnée de l'abscisse, sur un nombre impair de colonnes;
b) recherche de la droite la plus dense passant par le point précédent, par balayage du secteur [-45° ; +45°] et intervention de la fonction de pondération (Coef(d)) déjà mentionnée.

Ce que j'ai partiellement programmé fonctionne, mais il me faut reprendre l'orientation de l'algorithme.
Astuce éventuellement utile: afin de n'avoir à traiter qu'un seul type de fichier de dimensions déterminées, j'ai inséré l'image de la grille dans une autre image plus grande (600x500) - ici représentée en format réduit.

Le code précédent laissait un peu à désirer:

Les graphes semblent conduire à la valeur (X₀ = 7.5/2 = 3.75); l'ordonnée (Y₀) vaut toujours 73 .
Là réside le défaut de cette démarche exploratoire: il faut en fait obtenir pour (Y₀) une valeur moyenne (donc décimale) correspondant à une valeur entière de (X₀), arbitrairement choisie.

Je ne pourrai pas poursuivre avant quelques jours. Le sujet est intéressant, par les nouveaux algorithmes qu'il donne l'occasion de découvrir, et de mettre au point.

[wiwaxia]

J'ai regardé ce que donnait la sommation pondérée des indices de couleur à l'intérieur d'un cercle de rayon (Dmax), les coordonnées (x, y) du centre appartenant au domaine [Dmax ; La - Dmax]×[Dmax ; Ha - Dmax] .
Le coefficient de pondération dépend de la distance euclidienne (d) par la relation déjà indiquée:

Coef(d) = (1 - (d/Dmax)²)⁴ .

Deux représentations de la grille sont à ce stade envisageables, dont on peut saisir l'évolution au fur et à mesure que le rayon (Dmax) augmente:

a) une image en niveaux de gris, dont les lignes s'épaississent avant de devenir de plus en plus floues: résultat inexploitable en ce qui concerne la localisation des noeuds;

b) une image en fausses couleurs, qui présentera le contraste maximal à condition que soient connues les valeurs extrêmes de la somme (s) en chaque point (vert par exemple pour s = Smin, rouge pour s = Smax):

Les noeuds finissent bien par apparaître, mais avec une netteté très variable, et les éventuels défauts de l'image initiale dégradent la réponse attendue; le recours à une autre palette (FUNCTION F_Couleur_1(t: Reel): Pixel; ), dans laquelle les couleurs extrêmes sont plus accentuées, conduit à une meilleure localisation des noeuds, sans éliminer toutefois les défauts précédents.

Ci-dessous le coeur du programme Pascal calculant, à partir de la matrice du corps de l'image initiale (Matrice_0), celles des deux images successives (Matrice_1) d'abord en niveau de gris puis en fausses couleurs.
La première étape comporte la détermination des valeurs extrêmes (Smin, Smax) de la somme pondérée.

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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
 ... / ...
 TYPE Pixel = ARRAY[1..3] OF Byte;
      Tab_Pix = ARRAY[0..Dim_Max, 0..Dim_Max] OF Pixel;
 
 VAR Larg_Image, Haut_Image, T_Fichier, T_Image,
     Xmax, Xzero, Ymax, Yzero: Z_32;
     V_Fich_0, V_Fich_2, V_Fich_3: File OF Byte;
     Matrice_0, Matrice_1: Tab_Pix;
... / ...
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 P2 / Calcul de la seconde matrice
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 PROCEDURE ZeroM(Pfond: Pixel; VAR Ma: Tab_Pix);
   VAR i, j: Z_32;
   BEGIN
     FOR i:= 0 TO Dim_Max DO
       FOR j:= 0 TO Dim_Max DO Ma[i, j]:= Pfond
   END;
 
 PROCEDURE AffSminSmax(S_min, S_max: Reel);
   CONST C1 = 4; L1 = 44;
   BEGIN
     E(0015); Wt(C1, L1 - 1, 'Valeur minimale de la somme:    Smin = ');
     E(0012); Write(S_min:10:4);
     E(0015); Wt(C1, L1 + 1, 'Valeur maximale de la somme:    Smax = ');
     E(0012); Write(S_max:10:4)
   END;
 
 FUNCTION F_Couleur_1(t: Reel): Pixel;    // t varie sur [0 ; 1]
   CONST M = 255.499;
   VAR u, v, w: Reel; Px: Pixel;
   BEGIN
     w:= Sqr(1 - t);        u:= Sqr(1 - w);       Px[1]:= Round(M * u);
     w:= Sqr(t);            v:= Sqr(1 - w);       Px[2]:= Round(M * v);
     w:= (4 * t) * (1 - t); Px[3]:= Round(M * w); Result:= Px
   END;
 
 FUNCTION F_Couleur(t: Reel): Pixel;    // t varie sur [0 ; 1]
   CONST M = 255.499;
   VAR u, v, w: Reel; Px: Pixel;
   BEGIN
     u:= Sqr(t);            Px[1]:= Round(M * u);
     v:= Sqr(1 - t);        Px[2]:= Round(M * v);
     w:= (4 * t) * (1 - t); Px[3]:= Round(M * w); Result:= Px
   END;
 
 FUNCTION Coef(i, j, L2: Z_32): Reel;
   VAR s: Z_32; v, w: Reel;             // C = (1 - dý/Dý)^4
   BEGIN
     s:= Sqr(i); Inc(s, Sqr(j));
     IF (s>L2) THEN w:= 0
               ELSE BEGIN
                      v:= s / L2; w:= Sqr(1 - v)
                    END;
     Result:= Sqr(w)
   END;
 
 CONST D1 = 1;
       Pfond1: Pixel = (100, 100, 255);
       Pfond2: Pixel = (100, 255, 100);
 
 PROCEDURE Calc_Mat_Im2_B(La, Ha, X0, Y0: Z_32;
                          S_max, S_min: Reel; VAR Ma0, Ma1: Tab_Pix);
   CONST D2 = D1 * D1;
   VAR SIc, u, v, x, X1, X2, y, Y1, Y2: Z_32;
       g, Ks, s: Reel; Px: Pixel;
   BEGIN
     ZeroM(Pfond2, Matrice_1); Ks:= 1 / (S_max - S_min);
     X1:= D1;                  X2:= X0 - D1;
     Y1:= D1;                  Y2:= Y0 - D1;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         BEGIN
           s:= 0;
           FOR u:= -D1 TO D1 DO
             FOR v:= -D1 TO D1 DO BEGIN
                                    Px:= Ma0[x + u, y + v]; SIc:= Px[1];
                                    Inc(SIc, Px[2]);        Inc(SIc, Px[3]);
                                    g:= Coef(u, v, D2);     IncR(s, g * SIc)
                                  END;
           Ma1[x, y]:= F_Couleur(Ks * (s - S_min))
         END
   END;
 
 PROCEDURE Calc_Mat_Im2_A(La, Ha, X0, Y0: Z_32;
                          VAR S_max, S_min: Reel; VAR Ma0, Ma1: Tab_Pix);
   CONST D2 = D1 * D1;
   VAR k: Byte; SIc, u, v, x, X1, X2, y, Y1, Y2: Z_32;
       g, Ks, s, Max, Min: Reel; Px: Pixel;
   BEGIN
     ZeroM(Pfond1, Matrice_1); Min:= 1E9;    Max:= 0;
     X1:= D1;                  X2:= X0 - D1;
     Y1:= D1;                  Y2:= Y0 - D1; Ks:= 0.533 / D2;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         BEGIN
           s:= 0;
           FOR u:= -D1 TO D1 DO
             FOR v:= -D1 TO D1 DO BEGIN
                                    Px:= Ma0[x + u, y + v]; SIc:= Px[1];
                                    Inc(SIc, Px[2]);        Inc(SIc, Px[3]);
                                    g:= Coef(u, v, D2);     IncR(s, g * SIc)
                                  END;
           IF (Min>s) THEN Min:= s; IF (Max<s) THEN Max:= s;
           g:= Ks * s;              IF (g>255) THEN g:= 255;
           FOR k:= 1 TO 3 DO Px[k]:= Round(g);
           Ma1[x, y]:= Px
         END;
     S_max:= Max; S_min:= Min
   END;

Précisions à toutes fins utiles: Z_32 désigne le type LongInt, Reel le type Extended;
la procédure IncR(u, v) incrémente le réel (u) de la quantité réelle (v).

Finalement, le procédé ne permet pas tant de localiser les noeuds du réseau que les centres des mailles, ce qui revient strictement au même pour la résolution du problème en question.
Une version simplifiée de cet algorithme pourra intervenir dans la solution que j'ai évoquée, avec les mêmes éléments de programmation.

[wiwaxia]

Le procédé décrit conduit finalement au résultat attendu, malgré l'aspect initialement décevant des images obtenues; celles-ci ne sont que la représentation des variations d'une fonction de deux variables, dont il faut repérer les extremums dans un secteur donné.

Pour des raisons pratiques, l'image étudiée (158x143) a encore été incorporée à une autre de plus grandes dimensions (750x500); l'origine des coordonnées se situant au coin inférieur gauche, le domaine de la grille est défini par:

0 <= x <= Xmax ; 0 <= y <= Ymax

dont les bornes supérieures

Xmax = 157 ; Ymax = 142

sont rapidement déterminées par un algorithme simple:

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
 PROCEDURE Aff_XYmax;
   CONST C1 = 4; L1 = 39; o = 4;
   BEGIN
     E(0015); Wt(C1, L1 - 1, 'Limites de l''image initiale:');
     Wt(C1 + 8 , L1 + 1, 'Xmax = ');
     E(0012); Write(Xmax:o, '   '); E(0015); Write('Ymax = ');
     E(0012); Write(Ymax:o)
   END;
 
 FUNCTION Egalite_Pix(Pa, Pb: Pixel): Bool;
   VAR T1, T2, T3: Bool;
   BEGIN
     T1:= (Pa[1]=Pb[1]); T2:= (Pa[2]=Pb[2]);
     T3:= (Pa[3]=Pb[3]); Result:= (T1 AND (T2 AND T3))
   END;
 
 PROCEDURE Calc_XYmax(La, Ha: Z_32;VAR X_m, Y_m: Z_32; VAR Ma: Tab_Pix);
   VAR k: Z_32; PxC: Pixel;
   BEGIN
     PxC:= Ma[La - 1, Ha - 1];
     k:= 0; REPEAT Inc(k) UNTIL Egalite_Pix(Ma[k, 0], PxC); X_m:= k - 1;
     k:= 0; REPEAT Inc(k) UNTIL Egalite_Pix(Ma[0, k], PxC); Y_m:= k - 1;
     Aff_XYmax
   END;

Les dimensions de l'image de la grille correspondent à une distance approximative de 9 à 10 carreaux (Ncar, Ncar + 1); (D1) représentant par ailleurs la distance maximale à laquelle les pixels interviennent dans la sommation effectuée, un noeud du réseau se trouve nécessairement dans le domaine: [X1 ; X2]×[Y1 ; Y2]
dont les bornes sont définies par les relations:

X1= D1 + Round(Xm / Ncar) ; X2= X1 + Round(Xm / Ncar)
z = Ym DIV 2 ; Y1 = z - Round(z / Ncar); Y2 = z + Round(z / Ncar);

La localisation de ce point découle de la recherche du minimum de la somme pondérée dont il a déjà été question:

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
 CONST D1 = 7; D2 = D1 * D1; Ncar = 9;
       Pfond1: Pixel = (100, 100, 255);
       Pfond2: Pixel = (100, 255, 100);
... / ...
 PROCEDURE AffXYcen;
   CONST C1 = 4; L1 = 49; o = 4;
   BEGIN
     E(0015); Wt(C1, L1 - 1, 'Coordonn‚es du minimum rep‚r‚ (pour D1 = ');
     E(0010); Write(D1:2); E(0015); Write('):');
     Wt(C1 + 8 , L1 + 1, 'Xcen = ');
     E(0012); Write(Xcen:o, '   '); E(0015); Write('Ycen = ');
     E(0012); Write(Ycen:o)
   END;
 
 PROCEDURE Calc_XYcen(Xm, Ym: Z_32; VAR X_c, Y_c: Z_32; VAR Ma0: Tab_Pix);
   VAR h, SIc, u, v, x, X1, X2, Xc, y, Y1, Y2, Yc: Z_32;
       g, s, Smin: Reel; Px: Pixel;
   BEGIN
     h:= Round(Xm / Ncar); X1:= D1 + h; X2:= X1 + h;
     Y1:= Ym DIV 2;        Y2:= Y1;     h:= Round(Y1 / Ncar);
     Dec(Y1, h);           Inc(Y2, h);  Smin:= 1E10;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         BEGIN
           s:= 0;
           FOR u:= -D1 TO D1 DO
             FOR v:= -D1 TO D1 DO BEGIN
                                    Px:= Ma0[x + u, y + v]; SIc:= Px[1];
                                    Inc(SIc, Px[2]);        Inc(SIc, Px[3]);
                                    g:= Coef(u, v, D2);     IncR(s, g * SIc)
                                  END;
           IF (Smin>s) THEN BEGIN
                              Smin:= s; Xc:= x; Yc:= y
                            END
         END;
     X_c:= Xc; Y_c:= Yc; AffXYcen
   END;

Il s'agit du centre (C) qui joue par la suite un rôle déterminant. Pour le représenter sur les images déjà obtenues (en niveaux de gris ou en couleur), il suffit d'insérer dans le calcul de la matrice de pixels la procédure

Croix_C(1, Xcen, Ycen, Matrice_1) ou Croix_C(2, Xcen, Ycen, Matrice_1) .

En dehors de ce détail, l'algorithme ne fait que reprendre les anciennes instructions:

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
 
 PROCEDURE AffSminSmax(S_min, S_max: Reel);
   CONST C1 = 4; L1 = 44;
   BEGIN
     E(0015); Wt(C1, L1 - 1, 'Valeur minimale de la somme:      Smin = ');
     E(0012); Write(S_min:10:4);
     E(0015); Wt(C1, L1 + 1, 'Valeur maximale de la somme:      Smax = ');
     E(0012); Write(S_max:10:4)
   END;
 
 FUNCTION F_Couleur(t: Reel): Pixel;    // t varie sur [0 ; 1]
   CONST M = 255.499;
   VAR u, v, w: Reel; Px: Pixel;
   BEGIN
     u:= Sqr(t);            Px[1]:= Round(M * u);
     v:= Sqr(1 - t);        Px[2]:= Round(M * v);
     w:= (4 * t) * (1 - t); Px[3]:= Round(M * w); Result:= Px
   END;
 
 FUNCTION Coef(i, j, L2: Z_32): Reel;
   VAR s: Z_32; v, w: Reel;             // C = (1 - dý/Dý)^4
   BEGIN
     s:= Sqr(i); Inc(s, Sqr(j));
     IF (s>L2) THEN w:= 0
               ELSE BEGIN
                      v:= s / L2; w:= Sqr(1 - v)
                    END;
     Result:= Sqr(w)
   END;
 
 PROCEDURE Croix_C(c: Byte; Xc, Yc: Z_32; VAR Ma1: Tab_Pix);
   CONST Coul1: Pixel = (255, 000, 000);
         Coul2: Pixel = (255, 255, 000);
   VAR z:Z_32; Px: Pixel;
   BEGIN
     IF (c=1) THEN Px:= Coul1
              ELSE Px:= Coul2;
     FOR z:= -D1 TO D1 DO BEGIN
                            Ma1[Xc + z, Yc]:= Px;
                            Ma1[Xc, Yc + z]:= Px
                          END
   END;
 
 PROCEDURE Calc_Mat_Im1_B(La, Ha, X0, Y0: Z_32;
                          S_max, S_min: Reel; VAR Ma0, Ma1: Tab_Pix);
   VAR SIc, u, v, x, X1, X2, y, Y1, Y2: Z_32;
       g, Ks, s: Reel; Px: Pixel;
   BEGIN
     ZeroM(Pfond2, Matrice_1); Ks:= 1 / (S_max - S_min);
     X1:= D1;                  X2:= X0 - D1;
     Y1:= D1;                  Y2:= Y0 - D1;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         BEGIN
           s:= 0;
           FOR u:= -D1 TO D1 DO
             FOR v:= -D1 TO D1 DO BEGIN
                                    Px:= Ma0[x + u, y + v]; SIc:= Px[1];
                                    Inc(SIc, Px[2]);        Inc(SIc, Px[3]);
                                    g:= CoefVarR(u, v, D2); IncR(s, g * SIc)
                                  END;
           Ma1[x, y]:= F_Couleur(Ks * (s - S_min))
         END;
     Croix_C(2, Xcen, Ycen, Matrice_1)
   END;
 
 PROCEDURE Calc_Mat_Im1_A(La, Ha, X0, Y0: Z_32;
                          VAR S_max, S_min: Reel; VAR Ma0, Ma1: Tab_Pix);
   VAR k: Byte; SIc, u, v, x, X1, X2, y, Y1, Y2: Z_32;
       g, Ks, s, Max, Min: Reel; Px: Pixel;
   BEGIN
     ZeroM(Pfond1, Matrice_1); Min:= 1E9;    Max:= 0;
     X1:= D1;                  X2:= X0 - D1;
     Y1:= D1;                  Y2:= Y0 - D1; Ks:= 0.533 / D2;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         BEGIN
           s:= 0;
           FOR u:= -D1 TO D1 DO
             FOR v:= -D1 TO D1 DO BEGIN
                                    Px:= Ma0[x + u, y + v]; SIc:= Px[1];
                                    Inc(SIc, Px[2]);        Inc(SIc, Px[3]);
                                    g:= Coef(u, v, D2);     IncR(s, g * SIc)
                                  END;
           IF (Min>s) THEN Min:= s; IF (Max<s) THEN Max:= s;
           g:= Ks * s;              IF (g>255) THEN g:= 255;
           FOR k:= 1 TO 3 DO Px[k]:= Round(g);
           Ma1[x, y]:= Px
         END;
     Croix_C(1, Xcen, Ycen, Matrice_1);
     S_max:= Max; S_min:= Min
   END;

Voici les deux images obtenues, dans le cas où D1 = 7 :

Les coordonnées ne dépendent pas de la valeur de (D1), tant que celle-ci reste relativement faible (3 à 9):

[wiwaxia]

Après le hors d'oeuvre, le plat de résistance .

Soit deux demi-cercles centrés au point (C), dont le plus grand, tangent au bord le plus proche de l'image, admet pour rayon

R_ext = Min(Y_c, Y_max - Y_c) ;

Celui du demi-cercle intérieur (CM₀ = R_max) a sa valeur arbitrairement fixée à (R_ext - D₁),
et l'abscisse (X_C) du centre (C) est supérieure à la distance choisie (D₁).
On considère dans ce qui suit des segments dont l'inclinaison (t) est comprise entre (-90°) et (+90°).

Il résulte de ce qui précède que tout point (M) du segment (CM₀) est entouré de toutes parts des pixels de l'image, jusqu'à une distance au moins égale à (D₁).

On envisage en conséquence, pour un angle (t) donné, un certain nombre de points (M) (ici J_max = 200) régulièrement répartis sur (CM₀) et admettant pour coordonnées:

r = R_max * (j / J_max) ; x = X_C + r * Cos(t) ; y = Y_C + r * Sin(t) ;

il s'agit de points du domaine réel [0 ; Xm]×[0 ; Ym], dont les coordonnées sont généralement des nombres irrationnels, à la différence des pixels caractérisés par des coordonnées entières.
Parmi ces derniers, ceux dont la distance à (M) ne dépasse pas (D₁) se situent dans le domaine [X1 ; X2]×[Y1 ; Y2] défini par les relations:

X1 = Trunc(x) - D1; X2 = Trunc(x) + (D1 + 1);
Y1 = Trunc(y) - D1; Y2 = Trunc(y) + (D1 + 1);

d'où l'expression de la somme pondérée des indices de couleur, pour chacun des points (M), et la somme totale:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
         FOR j:= 1 TO Jmax DO          // Enumération des points successifs (M)
           BEGIN
             r:= Dr * j;
             x:= Som_2R(Xcen, r * Ct);    y:=  Som_2R(Ycen, r * St);
             Zm:= Trunc(x); X1:= Zm - D1; X2:= Zm + (D1 + 1);          // Calcul des limites du domaine
             Zm:= Trunc(y); Y1:= Zm - D1; Y2:= Zm + (D1 + 1);          
             FOR u:= X1 TO X2 DO
               FOR v:= Y1 TO Y2 DO BEGIN
                                     Px:= Ma0[u, v];  SIc:= Px[1];
                                     Inc(SIc, Px[2]); Inc(SIc, Px[3]);        // Somme des indices de couleur
                                     g:= CoefVarR(u - x, v - y, D2);          // Calcul du coefficient de pondération
                                     IncR(s, g * SIc)                         // Incrémentation de la somme (s)
                                   END
           END;

Deux différences, par rapport aux procédures précédentes:
a) l'intervention d'une fonction coefficient dépendant désormais de variables réelles (FUNCTION CoefVarR(Dx, Dy, R2: Reel): Reel;);
b) la consignation des (2*N_pas + 1) valeurs obtenues dans un tableau, afin de permettre la détermination indispensable des valeurs extrêmes rencontrées (S_min, S_max) lorsque l'inclinaison (t) varie de -90° à +90° :

Code :

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

1
2
3
4
5
6
7
8
9
10
   BEGIN
     IF ((2 * Yc)<Ym) THEN Rmax:= Yc ELSE Rmax:= Ym - Yc;
     Dec(Rmax, D1); Dr:= Rmax / Jmax;
     FOR i:= -Npas TO Npas DO
       BEGIN
         r:= i / Npas;      t:= H_Pi * r;
         SinCos(t, St, Ct); s:= 0;
         FOR j:= 1 TO Jmax DO
... / ...
         Liste[i]:= s                 END

Le calcul de la matrice de la dernière image n'intervient qu'ensuite - ci-dessous la partie du programme concernée:

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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
 CONST Npas = 360;
 
 TYPE LstR = ARRAY[-Npas..Npas] OF Reel;
 
 VAR Liste: LstR;
 
 PROCEDURE Aff_Min(J1, J2: Z_32);
   CONST C1 = 4;  C2 = C1 + 26; C3 = C1 + 32; C4 = C1 + 43;
         L1 = 54; v = 3;        w = 803;
   VAR Kt, T1, T2: Reel;
   BEGIN
     Kt:= 90 / Npas; T1:= Kt * J1; T2:= Kt * J2;
     E(0015);        Wt(C1, L1, 'Inclinaison des droites:');
     Wt(C2, L1 - 1, 'J1 =        T1 = ');
     Wt(C2, L1 + 1, 'J2 =        T2 = ');
     E(0012); We(C3, L1 - 1, J1, v); We(C3, L1 + 1, J2, v);
     E(0010); Wr(C4, L1 - 1, T1, w); Wr(C4, L1 + 1, T2, w)
   END;
 
 PROCEDURE Trace_Minima(Ha, Xm: Z_32; Lim: Reel; Li: LstR;
                        VAR J_1, J_2: Z_32; VAR Ma2: Tab_Pix);
   CONST P100: Pixel = (255, 0, 0);
   VAR j, J1, J2: Z_32; s, Sa, Sb: Reel;
   BEGIN
     J1:= 2 * Npas; J2:= 2 * Npas;
     FOR j:= (1 - Npas) TO (Npas - 1) DO
       BEGIN
         s:= Li[j];
         IF (s<Lim) THEN BEGIN
                            Sa:= Li[j - 1]; Sb:= Li[j + 1];
                            IF ((Sa>s) AND (s<Sb)) THEN
                              IF (J1>Npas) THEN J1:= j
                                           ELSE J2:= j
                         END
       END;
     FOR j:= 0 TO (Ha DIV 2) DO BEGIN
                                  Ma2[Xm + J1, j]:= P100;
                                  Ma2[Xm + J2, j]:= P100
                                END;
     J_1:= J1; J_2:= J2
   END;
 
 PROCEDURE Fond_Graphique(La, Ha, X_1, X_2, X_m, Y_1, Y_2: Z_32;
                          VAR Ma2: Tab_Pix);
   CONST m = 255;
         P000: Pixel = (0, 0, 0);
         P010: Pixel = (0, m, 0);
         P101: Pixel = (180, 000, 180);
         P111: Pixel = (m, m, m);
   VAR x, y: Z_32; TestX, TestY: Bool; Px: Pixel;
   BEGIN
     FOR x:= 0 TO (La - 1) DO
        FOR y:= 0 TO (Ha - 1) DO
          BEGIN
            TestX:= ((x<X_1) OR (x>X_2));
            TestY:= ((y<Y_1) OR (y>Y_2));
            IF (TestX OR TestY) THEN Px:= P101 ELSE Px:= P000;
            Ma2[x, y]:= Px
          END;
     FOR y:= 0 TO (Ha - 1) DO BEGIN
                                Ma2[X_1, y]:= P111; Ma2[X_2, y]:= P111;
                                Ma2[X_m, y]:= P111
                              END;
   END;
 
 PROCEDURE Calc_Mat_Im2(La, Ha: Z_32; Li: LstR; VAR Ma2: Tab_Pix);
   CONST P010: Pixel = (0, 255, 0);
   VAR h, Hmax, i, I1, I2, x, X1, X2, Xmed, y, Y1, Y2: Z_32;
       Ds, Kh, s, Slim, Smax, Smin: Reel; Px: Pixel;
   BEGIN
     Smax:= 0; Smin:= Li[1];
     FOR i:= -Npas TO Npas DO
       BEGIN
         s:= Li[i]; IF (Smax<s) THEN Smax:= s;
         IF (Smin>s) THEN Smin:= s
       END;
     Xmed:= La DIV 2;        X1:= Xmed - Npas;      X2:= Xmed + Npas;
     Hmax:= Round(0.9 * Ha); Y1:= Round(0.05 * Ha); Y2:= Y1 + Hmax;
     Kh:= Hmax / (Smax - Smin);
 
     Fond_Graphique(Larg_Image, Haut_Image, X1, X2, Xmed, Y1, Y2, Matrice_2);
 
     FOR i:= -Npas TO Npas DO
       BEGIN
         Ds:= Li[i] - Smin; h:= Round(Kh * Ds); x:= Xmed + i;
         FOR y:= Y1 TO (Y1 + h) DO Ma2[x, y]:= P010
       END;
     Slim:= Som_2R(0.8 * Smin, 0.2 * Smax);
     Trace_Minima(Haut_Image, Xmed, Slim, Liste,
                  Jmin_1, Jmin_2, Matrice_2);
     Aff_Min(Jmin_1, Jmin_2)
   END;
 
 FUNCTION CoefVarR(Dx, Dy, R2: Reel): Reel;
   VAR Dx2, Dy2, s, v, w: Reel;             // C = (1 - dý/Dý)^4
   BEGIN
     Dx2:= Sqr(Dx); Dy2:= Sqr(Dy); s:= Dx2 + Dy2;
     IF (s>R2) THEN w:= 0
               ELSE BEGIN
                      v:= s / R2; w:= Sqr(1 - v)
                    END;
     Result:= Sqr(w)
   END;
 
 PROCEDURE Calc_Liste(Ym, Yc: Z_32;
                      VAR L_: LstR; VAR Ma0: Tab_Pix);
   CONST Jmax = 200;
   VAR i, j, Rmax, SIc, u, v, X1, X2, Y1, Y2, Zm: Z_32;
       Ct, Dr, g, r, s, St, t, x, y: Reel; Px: Pixel;
   BEGIN
     IF ((2 * Yc)<Ym) THEN Rmax:= Yc ELSE Rmax:= Ym - Yc;
     Dec(Rmax, D1); Dr:= Rmax / Jmax;
     FOR i:= -Npas TO Npas DO
       BEGIN
         r:= i / Npas;      t:= H_Pi * r;
         SinCos(t, St, Ct); s:= 0;
         FOR j:= 1 TO Jmax DO
           BEGIN
             r:= Dr * j;
             x:= Som_2R(Xcen, r * Ct);    y:=  Som_2R(Ycen, r * St);
             Zm:= Trunc(x); X1:= Zm - D1; X2:= Zm + (D1 + 1);
             Zm:= Trunc(y); Y1:= Zm - D1; Y2:= Zm + (D1 + 1);
             FOR u:= X1 TO X2 DO
               FOR v:= Y1 TO Y2 DO BEGIN
                                     Px:= Ma0[u, v];  SIc:= Px[1];
                                     Inc(SIc, Px[2]); Inc(SIc, Px[3]);
                                     g:= CoefVarR(u - x, v - y, D2);
                                     IncR(s, g * SIc)
                                   END
           END;
         Liste[i]:= s
       END
   END;

Voici les graphes obtenus, correspondant aux distances maximales D₁ = 3, 5 et 7 ; les directions privilégiées de la grille s'y manifestent par deux pics orientés vers le bas:

On trouve (T₂ - T₁) = 90° , conformément au fait que la grille résulte de la superposition de deux faisceaux orthogonaux de droites parallèles; l'algorithme pourrait fonctionner avec tout autre réseau de points (monoclinique, hexagonal).

Les résultats sont légèrement affectés par la modification des paramètres (D₁, J_max); cette question est cependant indissociable de leur précision.
Il faudrait reprendre le tracé du graphe sur des intervalles plus étroits, centrés sur les valeurs déjà trouvées; l'amélioration de la précision n'est pas illimitée, car l'épaisseur des traits constitue un plancher en-dessous duquel il n'est pas possible de descendre.

[Guesset]

Bonjour,

Je trouve le problème et les réflexions très intéressants.

Juste pour apporter un autre axe de réflexion...

Comme nous avons un quadrillage a priori carré, nous avons 4 solutions tournées d'un quart de tour. Je supposerai que nous appliquions un principe d'économie en optant pour la solution d'angle le plus faible (les autres se déduisent facilement). Intéressons nous aux droites qui devraient être horizontales. Avec le principe d'économie, leur angle actuel doit être entre -45° et +45°.

A l'aide d'une verticale à gauche en Xo, cherchons un point optimal local Xo,Yo que nous considérerons a priori comme appartenant à une quasi-horizontale.

Si nous utilisons un Bresenham inversé, nous savons que pour X = X+1, l'optimal suivant ne peut être qu'en Y+1, Y, Y-1 (45° oblige). Prenons comme solution le Y qui est le plus optimal.

Cherchons à rejoindre le coté droit en suivant ce processus (c'est la balade dans le canyon).

Si cela échoue, c'est que l'hypothèse que le point Xo, Yo appartienne à une quasi horizontale est fausse. Dans ce cas (vite détecté) on prend comme nouveau point Xo, Yo, l'intersection suivante avec la verticale initiale. Un deuxième échec implique de décaler la verticale de telle sorte que (Xo+dX, Yo) soit hors du quadrillage (dans un carré) et de réitérer (cela signifie que l'angle cherché est très faible).

Si cela réussit, le point final Xf, Yf donne la pente cherchée : Xf-Xo, Yf-Yo.

Remarque : dès que |dY| est supérieur à une certaine valeur nous ne pouvons plus avoir de doute sur le quadrant de la droite et donc ne tester que deux des trois options pour Y.

J'arrête là, le chat a besoin du clavier pour sa sieste.

Salutations