Discussion :

[nant44]

Bonjour,
J'aimerai calculer "indirectement" la matrice de distribution gaussienne 2D a partir de la distribution gaussienne 1D. le mot "indirectement" veut dire : sans utiliser la formule "exp(x^2 + y^2)".
voiçi mon code si qqn peut trouver la solution.
D'avance merci

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
 
clear all
close all
clc
 
% generation d'une gaussienne 1D
N = 7;
nN = (N-1)/2;
sigm = 4;
xx = -nN:nN;
arg1 = -(xx.*xx)/(2*sigm*sigm);
h1 = exp(arg1);
figure;
plot(xx,h1);grid
hold on
 
% calcul de la matrice d'indexation 
[xin,yin]=meshgrid(-nN:nN);
zin=sqrt(xin.^2+yin.^2);
[inutile,inutile,zin(:)]=unique(zin(:));
nmax = max(max(zin));
 
%interpolation de la gaussienne 1D suivant le nombre "nmax"
% Note : cette interpolation sert a calculer les valeurs correspondants aux
% indices de la matrice "zin". C'est vraisemblablement l'origine du
% probleme, mais je n'arrive pas à voir comment calculer la valeur de la
% gaussienne pour tous les indices de "zin" !
hu = h1(end-nN:end);
nn = linspace(0,nN,nmax);
y = interp1(0:nN,hu,nn);
plot(nn,y,'r');
 
% remplissage de la matrice de distribution gaussienne
ZZ=y(zin);
figure;
mesh(xin,yin,ZZ);
hold on
 
% comparaison avec la méthode directe 
arg = -(xin.*xin + yin.*yin)/(2*sigm*sigm);
h = exp(arg);
figure
mesh(xin,yin,h)
hold on

[mr_samurai]

salut,

je ne comprend pas l'utilité de ton code.
Le but final est de tracer un gaussienne ? ou d'avoir une distribution gaussienne ?

Si tu veux une distribution gaussienne multi-dimensionnelle regarde plutot du coté de la fonction normrnd .

[nant44]

le but n'est pas de générer une gaussienne multi-dimensionnelles mais plutot une distribution multi-dimensionnelles quelconque . L'usage de la gaussienne ici sert tout simplement à vérifier les résultats trouvées.

[nant44]

j'ai réussi à débloquer partiellement le problème.
Comme on pourra le voir en compilant le code ci-desous que la surface calculée interpole bien la surface exacte sauf quand la distance (zin) dépasse le paramètre nN. C'est logique parce que on ne peut pas interpoler une courbe au delà de ses extrémités.
Donc, je serai reconnaissant si qqn peut me donner une piste ...

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
 
clear all
close all
clc
 
% calcul approché d'une gaussienne 2D
 
% 1) generation d'une gaussienne 1D
N = 9;
nN = (N-1)/2;
sigm = 4;
xx = -nN:nN;
arg1 = -(xx.*xx)/(2*sigm*sigm);
h1 = exp(arg1);
figure;
plot(xx,h1);grid
hold on
 
% 2) interpolation de la gaussienne 1D 
hu = h1(end-nN:end);
nn = linspace(0,nN,100*N);
y = interp1(0:nN,hu,nn);
plot(nn,y,'r');
 
[xin,yin]=meshgrid(-nN:nN);
zin=sqrt(xin.^2+yin.^2);
ZZ = zeros(size(zin));
[L,C]=size(zin);
for ii=1:L
    for jj=1:C
        zzin = abs(zin(ii,jj) - nn(1:end));
        ind = find(zzin == min(zzin));
        ind = ind(1);
        ZZ(ii,jj)=y(ind);
    end
end
figure
mesh(xin,yin,ZZ)    
hold on
 
% 3) comparaison avec la gaussienne 2D exacte 
arg = -(xin.*xin + yin.*yin)/(2*sigm*sigm);
h = exp(arg);
mesh(xin,yin,h)
 
norm(h-ZZ,'fro')

[Jerome Briot]

Je ne suis pas sûr de bien suivre ton raisonnement...

Ton problème vient du fait que les deux surfaces ne correspondent pas aux extrémités, c'est ça ?

Si c'est le cas, c'est tout à fait normale :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
>> [min(h1) max(h1)]
 
ans =
 
    0.6065    1.0000
 
>> [min(h(:)) max(h(:))]
 
ans =
 
    0.3679    1.0000

Ou visuellement :

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
clear all
close all
clc
 
% calcul approché d'une gaussienne 2D
 
% 1) generation d'une gaussienne 1D
N = 9;
nN = (N-1)/2;
sigm = 4;
xx = -nN:nN;
arg1 = -(xx.*xx)/(2*sigm*sigm);
h1 = exp(arg1);
figure;
plot3(xx,zeros(1,numel(xx)),h1,'b-s','linewidth',2);grid
hold on
 
% 2) interpolation de la gaussienne 1D 
hu = h1(end-nN:end);
nn = linspace(0,nN,100*N);
nnn = linspace(0,nN,numel(hu))
y = interp1(nnn,hu,nn)
plot3(nn,zeros(1,numel(nn)),y,'r-o','linewidth',3);
 
[xin,yin]=meshgrid(-nN:nN);
zin=sqrt(xin.^2+yin.^2);
ZZ = zeros(size(zin));
[L,C]=size(zin);
for ii=1:L
    for jj=1:C
        zzin = abs(zin(ii,jj) - nn(1:end));
        ind = find(zzin == min(zzin));
        ind = ind(1);
        ZZ(ii,jj)=y(ind);
    end
end
% figure
mesh(xin,yin,ZZ)    
hold on
 
% 3) comparaison avec la gaussienne 2D exacte 
arg = -(xin.*xin + yin.*yin)/(2*sigm*sigm);
h = exp(arg);
mesh(xin,yin,h)
 
norm(h-ZZ,'fro')
 
axis vis3d

Ou alors j'ai manqué quelque chose ?

[nant44]

Merci Dut pour ta réponse.

Ton problème vient du fait que les deux surfaces ne correspondent pas aux extrémités, c'est ça ?

Oui c'est ca le problème, mais y a t il une solution pour résoudre ce problème à ton avis ?

[mr_samurai]

salut,

je pense qu'il faut tronquer à sqrt(2) la taille du vecteur du départ pour ne pas avoir de NaN ... Problème classique de rotation de matrice.

Les modifications sont en %<Changement> et %</Changement> :

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
 
 
clear all
close all
clc
 
% calcul approché d'une gaussienne 2D
 
% 1) generation d'une gaussienne 1D
N = 15;
nN = (N-1)/2;
sigm = 4;
xx = -nN:nN;
arg1 = -(xx.*xx)/(2*sigm*sigm);
h1 = exp(arg1);
figure;
plot3(xx,zeros(1,numel(xx)),h1,'b-s','linewidth',2);grid
hold on
 
% 2) interpolation de la gaussienne 1D 
hu = h1(end-nN:end);
nn = linspace(0,nN,100*N);
nnn = linspace(0,nN,numel(hu));
y = interp1(nnn,hu,nn);
plot3(nn,zeros(1,numel(nn)),y,'g-o','linewidth',3);
 
%---- <Changement> ------------------
nN2 = round(nN/sqrt(2))-1;
[xin,yin]=meshgrid(-nN2:nN2);
%---- </Changement> ------------------
zin=sqrt(xin.^2+yin.^2);
ZZ = zeros(size(zin));
[L,C]=size(zin);
for ii=1:L
    for jj=1:C
        zzin = abs(zin(ii,jj) - nn(1:end));
        ind = find(zzin == min(zzin));
        ind = ind(1);
        ZZ(ii,jj)=y(ind);
    end
end
% figure
mesh(xin,yin,ZZ)    
hold on
 
% 3) comparaison avec la gaussienne 2D exacte 
arg = -(xin.*xin + yin.*yin)/(2*sigm*sigm);
h = exp(arg);
mesh(xin,yin,h)
 
norm(h-ZZ,'fro')
 
axis vis3d

L'erreur d'approximation dans ce cas : 0.0173

++ et bon chance

[nant44]

Merci mr_samurai pour ta solution, je vais analyser l'analyser pour les différentes distributions que j'utilise et je te tiens au courant.

[Jerome Briot]

Au passage, il faut la plupart du temps éviter de faire ceci :

Code :

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

idx = find(X == min(X))

Et préférer :

Code :

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

1
2
 
[pasbesoin,idx] = min(X);

Le code sera plus performant

Donc ton code devient :

Code :

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

1
2
3
4
5
6
7
for ii=1:L
    for jj=1:C
        zzin = abs(zin(ii,jj) - nn(1:end));
        [pasbesoin,ind] = min(zzin);
        ZZ(ii,jj)=y(ind);
    end
end

Note : on peut aussi faire ceci :

Code :

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

1
2
 
[idx,idx] = min(X);

Mais je trouve personnellement que cela nuit à la lisibilité du code (sans aucun avantage en terme de temps d'exécution)