Discussion :

[millie]

Merci. Avec le 40%, avec le portage que j'ai fait en C++, j'obtiens ça :



Un grand merci

Je ne vais pas poster maintenant le code source C++ car il s'intègre en fait à une bibliothèque que je suis en train de réaliser. Et mettre le code tout seul là serait inutile sans le reste de la bibliothèque.

Bon, pour le fun, je mets le code du fichier, mais il est assez sale (des tests de cohérence non effectués et des exceptions pas vérifiées), je voulais dans un premier temps faire marcher la chose :

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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
#include "InpaintOperator.hpp"
 
#include <cmath>
 
using namespace FImage;
 
static float absolue(float a);
 
inline float absolue(float a)
{
  if(a<0)
    return -a;
  else
    return a;
}
 
 
void  InpaintOperator::compute(Image & out, const Image & c,  Image& mask)
{
  int width = c.getWidth();
  int height = c.getHeight();
  static Image * pNewmask = NULL;
 
  if(out.getNumComponents() != c.getNumComponents())
   throw std::invalid_argument("InpaintOperator::compute");
 
  if(mask.getHeight() != c.getHeight() || mask.getWidth() != c.getWidth())
   throw std::invalid_argument("InpaintOperator::compue");
 
  out.resize(c.getWidth(), c.getHeight());
 
  if(pNewmask == NULL)
  {
    pNewmask = new Image(mask);
  }
  else
  {
    if(mask.getNumComponents() != pNewmask->getNumComponents())
    {
      delete pNewmask;
      pNewmask = new Image(mask);
    }
  }
 
  Image &newmask = *pNewmask;
 
  out = c;
  newmask = mask;
 
  for(int canal = 0; canal < c.getNumComponents(); canal++)
  {
 
    /*taille des vecteurs dx et dy*/
    int n = 8;
    int dx [] =
      {
        -1,0,1,1,1,0,-1,-1
      };
    int dy []  =
      {
        -1,-1,-1,0,1,1,1,0
      };
 
    std::vector<int> contourInt1;
    std::vector<int> contourInt2;
 
    while(true)
    {
 
      // front-line masked/unmasked
      for (int y = 0; y < height; y++)
      {
        for (int x = 0; x < width; x++)
        {
          if (mask.getPixel(x,y, canal)==0)
            continue;
          for (int i = 0; i < n; i++)
          {
            int xk = x + dx[i];
            int yk = y + dy[i];
            if (xk<0 || xk>=width)
              continue;
            if (yk<0 || yk>=height)
              continue;
            if (mask.getPixel(xk,yk, canal)>0)
              continue;
            contourInt1.push_back(x);
            contourInt2.push_back(y);
            break;
          }
        }
      }
 
      // exit when no front-line
      if (contourInt1.empty())
        break;
 
      // isophotes continuation
      for(int j=0;j<contourInt1.size();j++)
      {
        int x = contourInt1[j];
        int y = contourInt2[j];
 
        float value=0;
        float wsum=0;
        for (int i = 0; i < n; i++)
        {
          int xk = x + dx[i];
          int yk = y + dy[i];
          if (xk<0 || xk>=width)
            continue;
          if (yk<0 || yk>=height)
            continue;
          if (mask.getPixel(xk,yk, canal)>0)
            continue;
 
          // gradient
          float norme = 0;
          float angle = 0;
          InpaintOperator::gradient(c, xk, yk, canal, norme, angle);
 
          // gradient normal = isophote direction
          angle+=M_PI/2;
 
          // Weight of the propagation:
 
          // 1. dotproduct ( gradient normal . propagation vector )
          float pscal = cos(angle)*(-dx[i]) + (-sin(angle))*(-dy[i]);
          pscal/=sqrt(dx[i]*dx[i]+dy[i]*dy[i]);
 
          // 2. gradient magnitude  (O -> omnidirectionnal)
          float w = (norme)*absolue(pscal)+(1-norme)*1;
 
          value += w*c.getPixel(xk,yk, canal);
          wsum+=w;
        }
        if (wsum<=0)
          continue;
        value/=wsum;
 
        // set new value
        out.setPixel(x,y, canal, (int)value);
 
        // pixel becomes unmasked
        newmask.setPixel(x,y,canal, 0);
      }
 
 
      for(int j = 0; j<height; j++)
        for(int i = 0; i<width; i++)
          mask.setPixel(i,j, canal, newmask.getPixel(i,j, canal));
 
      /*vide la lsite*/
      contourInt1.clear();
      contourInt2.clear();
    }
  }
 
}
 
 
 
 
void InpaintOperator::gradient(const Image& c, int x, int y, int canal, float & norme, float & angle)
{
  int width = c.getWidth();
  int height = c.getHeight();
 
  float cst1 = (0.25*(2-sqrt(2.0)));
  float cst2 = (0.5f*(sqrt(2.0)-1));
 
  int px = x-1;
  int nx = x+1;
  int py = y-1;
  int ny = y+1;
  if (px<0)
    px=0;
  if (nx>=width)
    nx=width-1;
  if (py<0)
    py=0;
  if (ny>=height)
    ny=height-1;
 
  float Ipp=c.getPixel(px,py, canal);
  float Ipc=c.getPixel(px,y, canal) ;
  float Ipn=c.getPixel(px,ny, canal);
  float Icp=c.getPixel(x,py, canal);
  float Icn=c.getPixel(x,ny, canal);
  float Inp=c.getPixel(nx,py, canal);
  float Inc=c.getPixel(nx,y, canal) ;
  float Inn=c.getPixel(nx,ny, canal);
 
  float IppInn = cst1*(Inn-Ipp);
  float IpnInp = cst1*(Ipn-Inp);
  float gradx = (IppInn-IpnInp-cst2*Ipc+cst2*Inc);
  float grady = (IppInn+IpnInp-cst2*Icp+cst2*Icn);
 
  norme = sqrt( gradx*gradx + grady*grady );
 
  angle = 0;
  if (norme>0)
  {
    angle = acos(gradx/norme);
    if (grady>0)
      angle = 2*M_PI - angle;
  }
 
  norme/=255;
 
}

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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#include <cstdio>
#include <cmath>
#include <iostream>
 
#include "InpaintNoiseFilter.hpp"
 
#include "InpaintOperator.hpp"
 
using namespace FImage;
 
static float absolue(float a);
 
inline float absolue(float a)
{
  if(a<0)
    return -a;
  else
    return a;
}
 
 
static inline float max(float a, float b)
{
  if(a<b)
    return b;
  else
    return a;
}
 
static inline int sup(const void *  a,  const void *b)
{
  return  ( *(int*)a - *(int*)b );
 
}
 
void InpaintNoiseFilter::computeRoad(Image & mask, Image & c, int aperture, float coef)
{
  static float distance[256];
  static DataBuffer<int> dist(1);
 
  dist.resize((2*aperture+1)*(2*aperture+1));
 
  int width = c.getWidth();
  int height = c.getHeight();
 
 
  distance[0]=0;
  for(int i=1;i<256;i++)
  {
    float x = (float)i/255;
    float logr = log(x)/log(2);
    distance[i]=1+max(logr,-5)/5;
  }
 
 
  for(int canal=0; canal< c.getNumComponents(); canal++)
  {
 
    for (int y=0; y<height; y++)
    {
      // std::cerr<<"I";
      for (int x=0; x<width; x++)
      {
 
        // compute distances between pixels
        int average=0;
        int count=0;
        for(int dy=-aperture;dy<=aperture;dy++)
        {
          for(int dx=-aperture;dx<=aperture;dx++)
          {
            if ((x+dx)<0)
              continue;
            if ((x+dx)>=width)
              continue;
            if ((y+dy)<0)
              continue;
            if ((y+dy)>=height)
              continue;
 
            if (c.getPixel(x+dx,y+dy, canal)<0)
              continue;
 
            int i = (int) absolue(c.getPixel(x+dx,y+dy, canal)-c.getPixel(x,y, canal));
            dist[count]=(int)(255.0*distance[i]);
            count++;
 
            average+= (int) c.getPixel(x+dx,y+dy, canal);
          }
        }
        average/=count;
 
        // Noise estimation (ROAD at rank nmb_pixels/2)
        qsort(dist.getDataBuffer(), count, sizeof(int), sup);
 
        int road = 0;
        int rank=count/2;
        for(int i=0;i<rank;i++)
          road+=dist[i];
 
        // compute threshold (= variance)
        int etype=0;
        count=0;
        for(int dy=-aperture;dy<=aperture;dy++)
        {
          for(int dx=-aperture;dx<=aperture;dx++)
          {
            if ((x+dx)<0)
              continue;
            if ((x+dx)>=width)
              continue;
            if ((y+dy)<0)
              continue;
            if ((y+dy)>=height)
              continue;
 
            float e = c.getPixel(x+dx,y+dy, canal)-average;
            etype+= (int) (e*e);
            count++;
          }
        }
        etype=(int) sqrt(etype/count);
        int threshold = (int)(rank*etype);
 
        // pixel exceed threshold -> noise
 
        if ( (road>=(threshold*coef)) && (etype>0) )
        {
          mask.setPixel(x,y,canal, 255);
        }
        else
        {
          mask.setPixel(x,y,canal, 0);
        }
 
 
 
      }
    }
 
 
 
  }
}
 
void InpaintNoiseFilter::compute(Image & out, Image & c, int aperture, float coef, int itermax)
{
  Image mask(c.getWidth(), c.getHeight(), c.getNumComponents());
 
  mask.resize(c.getWidth(),c.getHeight());
  Image work(c.getNumComponents());
 
  work = c;
 
  for(int i=0;i<itermax;i++)
  {
    computeRoad(mask, c,aperture, coef);
 
    InpaintOperator::compute(out, work,mask);
    work = out;
  }
 
 
}

[pseudocode]

EDIT 2 : Non, je fus mauvaise langue, le filtre marche vraiment bien

Je viens de capter d'ou venait le probleme. Tu as pris l'image de mon Post pour faire le test. Sauf que j'avais reduit l'image de moitié, par interpolation bicubique

==> un superbe lissage du "bruit"

==> l'estimateur de ROAD a considéré que pratiquement toute l'image de gauche etait du bruit

==> l'inpaint a recréé une image presque de zero. D'ou la jolie gouache que tu obtiens.


EDIT: c'est fou comme le C ressemble au Java... A moins que ce ne soit l'inverse ... (Troll inside)

EDIT2: Pourquoi tu as des pixels blancs avec Lena40% ? Moi j'en ai aucun.

[millie]

Je ne sais pas. Peut être parce que je travaille avec des float et que la précision est différente en C++ ? (ou j'ai peut être supprimé une ligne de ton code sans faire attention, en me disant qu'elle servait à rien )

Oui, le code est ressemblant car je me suis inspiré de JAI pour faire la conception de ma biblio (les méthodes étaient au début en français, et je me suis dit que j'allais tout remettre en anglais, et j'ai pris le nom des méthodes java).

[PRomu@ld]

A noter que cet algo ne marche bien que pour les images avec un bruit Salt & Pepper (points noirs et blancs). Par exemple, pour ameliorer des images recues par fax ou mal scannées (c'est mon cas).

On ne parle donc bien que de bruits impulsionnels, dans ce cas, un filtre médian n'est-il pas plus efficace (en terme de performance) et acceptable (en terme de "qualité").

Pour les bruits gaussiens ça donne quoi ? (des exemples, des analyses, ...)

[millie]

Pour les bruits gaussiens ça donne quoi ? (des exemples, des analyses, ...)

J'essayerai ce soir (si j'ai le temps )

On ne parle donc bien que de bruits impulsionnels, dans ce cas, un filtre médian n'est-il pas plus efficace (en terme de performance) et acceptable (en terme de "qualité").

Le filtre médian donne ça :

[millie]

Avec 20% de bruit gaussien, j'ai ça :



à



(avec les mêmes paramètres)

Avec la méthode des hypersurfaces, j'ai seulement réussi à faire ça

[PRomu@ld]

Pour le filtre médian tu utilises quoi ? 3*3 ?

Mais finalement le bruit après filtrage est trop important par rapport à l'autre méthode, c'était donc une mauvaise idée , désolé.

[millie]

Pour le filtre médian tu utilises quoi ? 3*3 ?

Mais finalement le bruit après filtrage est trop important par rapport à l'autre méthode, c'était donc une mauvaise idée , désolé.

En fait, j'ai bouletisé. J'ai appliqué un filtre moyenneur, et pas médian

[PRomu@ld]

C'est pas grave, j'aurai du aussi m'en rendre compte, le résultat est pourtant typique d'un moyenneur sur un bruit impulsionnel.

Mais comme c'est toi l'expert en TDI je m'étais dit que je pouvait te faire confiance

[millie]

Alors, en utilisant le filtre median de gimp (désolé, je n'ai pas encore codé cette partie dans ma bibliothèque (et en plus je suis pas chez moi) ) :




Il faut avouer que la qualité de l'autre est mieux mais ça reste pas mal

[PRomu@ld]

Oui, c'est moins joli, mais ça reste relativement acceptable (enfin tout dépend du traitement qu'on désire faire derrière),

Je doute qu'avec un réhausseur (laplacien + 1 au centre) après le médian on puisse approcher une qualité semblable à celle de l'autre filtre.

[pseudocode]

On ne parle donc bien que de bruits impulsionnels, dans ce cas, un filtre médian n'est-il pas plus efficace (en terme de performance) et acceptable (en terme de "qualité").

Pour les bruits gaussiens ça donne quoi ? (des exemples, des analyses, ...)

C'est sur que le filtre médian est plus rapide et plus robuste que celui-la. J'ai poster ce filtre sans prétentions. Je n'ai pas l'intention de revolutionner le domaine du filtrage .

Cette technique ne marche ques pour les bruits impulsionnels. En effet, l'idée est de considerer que chaque pixel de l'image est soit pur, soit bruité: tout l'un, ou tout l'autre.

Sur une image avec des bruits gaussien, chaque pixel est altéré. On ne peut donc pas utiliser les estimateurs pour séparer les pixel purs et les pixels bruités. Il restera toujours des pixels bruités. La restauration par inpainting ne marchera pas (les isophotes sont rompus).

[PRomu@ld]

C'est sur que le filtre médian est plus rapide et plus robuste que celui-la. J'ai poster ce filtre sans prétentions. Je n'ai pas l'intention de revolutionner le domaine du filtrage

Je ne remettais pas en cause ton filtre, je voulais seulement avoir une comparaison par rapport à un filtre médian (puisqu'ils permettent d'éliminer le même type de bruit), c'est donc pour peser le pour et le contre de chacun de filtre afin d'aider certains si jamais ils ont besoin d'appliquer un filtre.

Force est de constater que ton filtre donne de plus "joli" résultats (en tout cas sur lena), il faudrait regarder sur d'autres types d'images plus particulières pour vérifier les résultats, mais en tout cas, tout ceci est une bonne idée.

[millie]

Code :

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

1
2
3
4
5
[
				newmask.set(x,y,0);
			}
 
			mask=newmask;

Tu es sûr de ton mask = newmask (a priori, il faudrait faire une recopie) ?

En tout cas, j'ai de gros problèmes pour faire passer le code avec des optimisations. C'est une erreur très bizarre.

[pseudocode]

Tu es sûr de ton mask = newmask (a priori, il faudrait faire une recopie) ?

En tout cas, j'ai de gros problèmes pour faire passer le code avec des optimisations. C'est une erreur très bizarre.

Non, je ne suis sûr de rien .

Hum... non, c'est bon. En fait, c'est surtout que le code est sale. Normalement je devrais utiliser tout le temps "newmask" dans la boucle "while(true)", et donc je pourrais faire sauter la ligne "mask=newmask".

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
/**
 * Fast (not accurate) Inpainting 
 * 
 * @param c Input data (luminosity/gray-level)
 * @param mask masked pixel to evaluate
 * @return Output data (luminosity/gray-level)
 */
private ByteProcessor inpaint(ByteProcessor c, ByteProcessor mask) {
	int width = c.getWidth();
	int height = c.getHeight();
 
	ByteProcessor c2 = new ByteProcessor(width,height);
	for (int y=0; y<height; y++)
		for (int x=0; x<width; x++)
			c2.set(x,y,c.get(x,y));
 
	ByteProcessor newmask = new ByteProcessor(width,height);
	for (int y=0; y<height; y++)
		for (int x=0; x<width; x++)
			newmask.set(x,y,mask.get(x,y));
 
	int n = 8;
	int[] dx = new int[] {-1,0,1,1,1,0,-1,-1};
	int[] dy = new int[] {-1,-1,-1,0,1,1,1,0};
 
	while(true) {
 
		// front-line masked/unmasked
		List contourInt = new ArrayList();
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				if (newmask.get(x,y)==0) continue;
				for (int i = 0; i < n; i++) {
					int xk = x + dx[i];
					int yk = y + dy[i];
					if (xk<0 || xk>=width)  continue;
					if (yk<0 || yk>=height)  continue;
					if (newmask.get(xk,yk)>0) continue;
					contourInt.add( new int[] {x,y} );
					break;
				}
			}
		}
 
		// exit when no front-line
		if (contourInt.isEmpty()) break;
 
		// isophotes continuation
		for(int j=0;j<contourInt.size();j++) {
			int[] pixel = (int[])contourInt.get(j);
			int x =  pixel[0];
			int y =  pixel[1];
 
			double value=0;
			double wsum=0;
			for (int i = 0; i < n; i++) {
				int xk = x + dx[i];
				int yk = y + dy[i];
				if (xk<0 || xk>=width)  continue;
				if (yk<0 || yk>=height)  continue;
				if (newmask.get(xk,yk)>0) continue;
 
				// gradient
				double[] grad = gradient(c,xk,yk);
				double norme = grad[0];
				double angle = grad[1];
 
				// gradient normal = isophote direction
				angle+=Math.PI/2;
 
				// Weight of the propagation:
 
				// 1. dotproduct ( gradient normal . propagation vector )
				double pscal = Math.cos(angle)*(-dx[i]) + (-Math.sin(angle))*(-dy[i]);
				pscal/=Math.sqrt(dx[i]*dx[i]+dy[i]*dy[i]);
 
				// 2. gradient magnitude  (O -> omnidirectionnal)
				double w = (norme)*Math.abs(pscal)+(1-norme)*1;
 
				value += w*c.get(xk,yk);
				wsum+=w;
			}
			if (wsum<=0) continue;
			value/=wsum;
 
			// set new value
			c2.set(x,y,(int)value);
 
			// pixel becomes unmasked
			newmask.set(x,y,0);
		}
	}
 
	return c2;
}

[pseudocode]

Avec la méthode des hypersurfaces, j'ai seulement réussi à faire ça

Dans ces cas la, je prefere le filtre anisotropique (laplace-beltrami):

[AIDA DIANE]

bonjour, je suis une novice en programmation. je dois réaliser une cartographie d'image IRM avec image j. je souhaite utiliser UnNoise pour débruiter mes images. Cependant, je ne comprends pas quelle est la différence entre faire le filtrage avec ROAD et Variance et comment se fait le choix des valeurs de window aperture, variance factor et iterations. merci beaucoup

[pseudocode]

Ce filtre a été conçu pour retirer du bruit impulsionnel. Je doute que cela soit le cas pour des images IRM.

Ouvrez une discussion dans le forum "Traitement d'images" et postez une de vos images afin que nous puissions vous aider.