1. #1
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    10 027
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 44
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 10 027
    Points : 15 578
    Points
    15 578

    Par défaut [Image] Filtre de diffusion Laplace-Beltrami pour ImageJ

    Puisque j'en ai parlé précédement, je le poste.

    Laplace-Beltrami operator (http://en.wikipedia.org/wiki/Laplacian)

    Anisotropic Feature-Preserving Denoising of Height Fields and Bivariate Data (www.multires.caltech.edu/pubs/gi2000.pdf)


    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
    212
    213
    214
    215
    216
    217
    218
    219
    220
    221
    222
    223
    224
    225
    226
    227
    228
    229
    230
    231
    232
    233
    234
    235
    236
    237
    238
    239
    240
    241
    242
    243
    244
    245
    246
    247
    248
    249
    250
    251
    252
    253
    254
    255
    256
    257
    258
    259
    260
    261
    262
    263
    264
    265
    266
    267
    268
    269
    270
    271
    272
    273
    274
    275
    276
    277
    278
    279
    280
    281
    282
    283
    284
    285
    286
    287
    288
    289
    290
    291
    292
    293
    294
    295
    296
    297
    298
    299
    300
    301
    302
    303
    304
    import ij.IJ;
    import ij.ImagePlus;
    import ij.gui.GenericDialog;
    import ij.plugin.filter.PlugInFilter;
    import ij.process.ByteProcessor;
    import ij.process.ColorProcessor;
    import ij.process.ImageProcessor;
    
    import java.awt.Color;
    
    /**
     * @author Xavier Philippeau
     * 
     * LaplaceBeltrami Plugin. 
     *
     * Anisotropic Feature-Preserving Denoising using the Laplace-Beltrami operator
     */
    public class LaplaceBeltrami_ implements PlugInFilter {
    	
    	private ImagePlus imgOrig = null;
    	private int itermax = 0;
    	private boolean autostop=false;
    	private double step=0;
    	
    	/**
    	* This method gets a reference to the image to be locally ranked
    	* and returns the filter capabilities.
    	* @param	arg	calls for the "about" information.
    	* @param	imp	This is the image to be processed
    	*
    	* @return	Returns DOES_8G + DOES_RGB or DONE.
    	*/
    	public int setup(String arg, ImagePlus imp) {
    		
    		// about...
    		if (arg.equals("about")) {
    			showAbout(); 
    			return DONE;
    		}
    
    		// else...
    		if (imp==null) return DONE;
    
    		//original image
    		this.imgOrig = imp ;
    		
    		// Configuration dialog.
    		GenericDialog gd = new GenericDialog("Parameters");
    		gd.addNumericField("Integration step",0.1,1);
    		gd.addNumericField("Iterations (-1: auto)",-1,0);
    		
    		while(true) {
    			gd.showDialog();
    			if ( gd.wasCanceled() )	return DONE;
    		
    			
    			this.step = (double) gd.getNextNumber();
    			this.itermax = (int) gd.getNextNumber();
    			
    			if (this.step<=0) continue;
    			if (this.itermax<-1) continue;
    			break;
    		}
    		gd.dispose();
    		
    		if (this.itermax==-1) {
    			this.autostop = true;
    			this.itermax=Integer.MAX_VALUE;
    		}
    		
    		return DOES_8G+DOES_RGB;
    	}
    
    	/** 
    	 * Filters use this method to process the image. If the
    	 * SUPPORTS_STACKS flag was set, it is called for each slice in
    	 * a stack. ImageJ will lock the image before calling
    	 * this method and unlock it when the filter is finished. 
    	 */
    	public void run(ImageProcessor ip) {
    		// project image to HSV color space
    		ByteProcessor work = null;
    		if (this.imgOrig.getBitDepth()==24) work=transformRGB(this.imgOrig);
    		if (this.imgOrig.getBitDepth()==8) work=transformGray(this.imgOrig);
    		
    		// show new image		
    		ImagePlus newImg = new ImagePlus(this.imgOrig.getTitle()+"_LaplaceBeltrami", this.imgOrig.getProcessor());
    		newImg.show();
    
    		// filter
    		ByteProcessor previous = work;
    		int previousdiff = Integer.MAX_VALUE;
    		for(int i=0;i<this.itermax;i++) {
    			IJ.showProgress((float)i/this.itermax);
    			IJ.showStatus("Iteration: "+i);
    			
    			// compute one step
    			work = LaplaceBeltramiFilter(work,this.step);
    
    			// project image to RGB color space
    			ImageProcessor newIp = null;
    			if (this.imgOrig.getBitDepth()==24) newIp=invTransformRGB(this.imgOrig,work);
    			if (this.imgOrig.getBitDepth()==8) newIp=invTransformGray(this.imgOrig,work);
    			newImg.setProcessor(null,newIp);
    			newImg.updateAndDraw();
    
    			// auto exit at idempotence
    			if (this.autostop) {
    				int diff=compare(work,previous);
    		    	if (diff<=0) break;
    		    	int rate = (previousdiff-diff);
    		    	if (rate<=0) {
    		    		this.itermax = this.itermax+1;
    		    	} else {
    		    		this.itermax = i+2+diff/rate;
    		    	}
    		    	previousdiff=diff;
    		    	previous=work;
    		    }
    		}
    		IJ.showStatus("Done");
    		IJ.showProgress(0);
    	}
    
    	// RGB -> HSV
    	private ByteProcessor transformRGB(ImagePlus img) {
    		ByteProcessor bp = new ByteProcessor(img.getWidth(),img.getHeight());
    		for (int y = 0; y < img.getHeight(); y++) {
    			for (int x = 0; x < img.getWidth(); x++) {
    				int[] lrgb = img.getPixel(x,y);
    				float[] hsb = Color.RGBtoHSB(lrgb[0],lrgb[1],lrgb[2],null);
    				bp.set(x,y,(int)(hsb[2]*255));
    			}
    		}
    		return bp;
    	}
    	
    	// GRAY -> HSV
    	private ByteProcessor transformGray(ImagePlus img) {
    		ByteProcessor bp = new ByteProcessor(img.getWidth(),img.getHeight());
    		for (int y = 0; y < img.getHeight(); y++) {
    			for (int x = 0; x < img.getWidth(); x++) {
    				int[] lrgb = img.getPixel(x,y);
    				bp.set(x,y,lrgb[0]);
    			}
    		}
    		return bp;
    	}
    	
    	// HSV -> RGB
    	private ImageProcessor invTransformRGB(ImagePlus img, ByteProcessor bp) {
    		ImageProcessor out = new ColorProcessor(img.getWidth(),img.getHeight());
    		for (int y = 0; y < img.getHeight(); y++) {
    			for (int x = 0; x < img.getWidth(); x++) {
    				int[] lrgb = img.getPixel(x,y);
    				float[] hsb = Color.RGBtoHSB(lrgb[0],lrgb[1],lrgb[2],null);
    				int rgb = Color.HSBtoRGB(hsb[0],hsb[1],bp.get(x,y)/255f);
    				out.set(x,y,rgb);
    			}
    		}
    		return out;
    	}
    
    	// HSV -> GRAY
    	private ImageProcessor invTransformGray(ImagePlus img, ByteProcessor bp) {
    		ImageProcessor out = new ByteProcessor(img.getWidth(),img.getHeight());
    		for (int y = 0; y < img.getHeight(); y++) {
    			for (int x = 0; x < img.getWidth(); x++) {
    				out.set(x,y,bp.get(x,y));
    			}
    		}
    		return out;
    	}
    
    	// About...
    	private void showAbout() {
    		IJ.showMessage("Laplace-Beltrami...","Laplace-Beltrami Filter by Pseudocode");
    	}
    
    	// Compare two images
    	private int compare(ByteProcessor c,ByteProcessor c1) {
    		int diff=0;
    		for (int y=0;y<c.getHeight();y++) {
    			for (int x=0;x<c.getWidth();x++) {
    				if (Math.abs(c.get(x,y)-c1.get(x,y))>0) diff++;
    			}
    		}
    		return diff;
    	}
    	
    	// compute cotangent of the 2 vectors u and v
    	private double cotangent(int ux,int uy,int uz, int vx,int vy,int vz) {
    		double uv  = ux*vx + uy*vy + uz*vz;
    		
    		int vpx = uz*vy - uy*vz;
    		int vpy = ux*vz - uz*vx;
    		int vpz = ux*vy - uy*vx;
    		double vp = Math.sqrt( vpx*vpx + vpy*vpy + vpz*vpz );
    		if (vp==0) return Double.MAX_VALUE;
    			
    		return Math.abs(uv)/vp;
    	}
    
    	// compute area of a 3D triangle
    	private double area(int xa,int ya,int va,int xb,int yb,int vb,int xc,int yc,int vc) {
    		int ux = xa-xb;
    		int uy = ya-yb;
    		int uz = va-vb;
    		
    		int vx = xc-xb;
    		int vy = yc-yb;
    		int vz = vc-vb;
    
    		int vpx = uz*vy - uy*vz;
    		int vpy = ux*vz - uz*vx;
    		int vpz = ux*vy - uy*vx;
    		
    		double area = Math.sqrt(vpx*vpx + vpy*vpy + vpz*vpz)/2;
    		return area;
    	}
    		
    	// Laplace-Beltrami operator (divergence of the gradient) 
    	public ByteProcessor LaplaceBeltramiFilter(ByteProcessor c, double dt) {
    		int width = c.getWidth();
    		int height = c.getHeight();
    
    		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};
    		
    		ByteProcessor c2 = new ByteProcessor(width,height);
    
    		// For each pixel of the image
    		for (int y=0; y<height; y++) {
    			for (int x=0; x<width; x++) {
    				
    				int z = c.get(x, y);
    				c2.set(x, y, z);
    				
    				// Laplace-Beltrami Operator
    				// -------------------------
    				//
    				// L[i] = 1/2A * Sum[ ( cotan(Alphaij) + cotan(Betaij) ) * ( V(j) - V(i) )  ]
    				//               j in N(i)
    				//
    				// V(x): Value (intensity,height) of pixel x
    				// N(x) : neighboors of pixel x
    				// A : Global Area (sum) of surrounding triangles
    				// Alphaij, Betaij : Opposite angles of triangles containing the edge (i,j)
    				//
    				double laplace=0; double area=0;
    				for(int i=0;i<n;i++) {
    					
    					// previous vertex
    					int xp = x+dx[(i+n-1) % n]; 
    					if (xp<0 || xp>=width) continue;
    					int yp = y+dy[(i+n-1) % n]; 
    					if (yp<0 || yp>=height) continue;
    					int zp = c.get(xp, yp);
    
    					// actual vertex
    					int xk = x+dx[i]; 
    					if (xk<0 || xk>=width) continue;
    					int yk = y+dy[i];
    					if (yk<0 || yk>=height) continue;
    					int zk = c.get(xk, yk);
    					
    					// next vertex
    					int xn = x+dx[(i+1) % n]; 
    					if (xn<0 || xn>=width) continue;
    					int yn = y+dy[(i+1) % n];
    					if (yn<0 || yn>=height) continue;
    					int zn = c.get(xn, yn);
    
    					// cotangent of opposite angles of the 2 triangles
    					double cotgt_alpha = cotangent(x-xp,y-yp,z-zp, xk-xp,yk-yp,zk-zp);
    					double cotgt_beta  = cotangent(x-xn,y-yn,z-zn, xk-xn,yk-yn,zk-zn);
    					
    					// area of the 2 triangles
    					double area_alpha = area(x,y,z, xp,yp,zp, xk,yk,zk);
    					double area_beta = area(x,y,z, xn,yn,zn, xk,yk,zk);
    					
    					// addd to iteration variable
    					laplace += (cotgt_alpha+cotgt_beta)*(zk-z);
    					area += area_alpha+area_beta;
    				}
    				laplace/=2*area;
    
    				// Compute new value:
    				//
    				// dI/dt ~= ( I(t+dt) - I(t) )  / dt
    				// ==> I(t+dt) ~= I(t) + dt * (dI/dt) ~= I(t) + dt * L[I(t)]
    				
    				int v=(int)Math.rint( c.get(x, y) + dt*laplace );
    				if (v<0) v=0;
    				if (v>255) v=255;
    				c2.set(x, y, v);
    			}
    		}
    		
    		return c2;
    	}
    }
    PRomu@ld : sources
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  2. #2
    Rédacteur

    Avatar de millie
    Profil pro
    Inscrit en
    juin 2006
    Messages
    6 937
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : juin 2006
    Messages : 6 937
    Points : 9 555
    Points
    9 555

    Par défaut

    Je l'ai implémenté (sauf que je fais une itération bourrine sans recalculer le nombre d'itération).

    Sauf que j'ai des résultats pas tip top, j'avais essayé avec dt = 1 et le nombre d'itération aux alentours de 20. Faut il beaucoup plus ?
    Je ne répondrai à aucune question technique en privé

  3. #3
    Rédacteur

    Avatar de millie
    Profil pro
    Inscrit en
    juin 2006
    Messages
    6 937
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : juin 2006
    Messages : 6 937
    Points : 9 555
    Points
    9 555

    Par défaut

    Non, j'avais fait n'importe quoi (non, pas tant que ça, mais j'avais oublié une ligne à un endroit).

    Ca marche nickel
    Je ne répondrai à aucune question technique en privé

  4. #4
    Rédacteur
    Avatar de pseudocode
    Homme Profil pro
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    10 027
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Âge : 44
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 10 027
    Points : 15 578
    Points
    15 578

    Par défaut

    Je suppose que la version C++/MillieLib booste plus que la version Java/ImageJ ?
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

Discussions similaires

  1. [image] Filtre Squelette pour ImageJ
    Par pseudocode dans le forum Contribuez
    Réponses: 47
    Dernier message: 12/08/2013, 10h41
  2. [Image] Filtre UnNoise pour ImageJ
    Par pseudocode dans le forum Contribuez
    Réponses: 37
    Dernier message: 07/03/2008, 16h23
  3. [Image] Filtre de Canny pour ImageJ
    Par pseudocode dans le forum Contribuez
    Réponses: 18
    Dernier message: 13/09/2007, 19h01
  4. [Image] Filtre UnNoise pour ImageJ
    Par pseudocode dans le forum Général Algorithmique
    Réponses: 10
    Dernier message: 03/04/2007, 23h38
  5. [Image] Filtre de Laplace
    Par scifire dans le forum 2D
    Réponses: 5
    Dernier message: 21/10/2005, 17h15

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo