Image Edge Enhancing Coherence Filter Toolbox: u=diffusion_scheme_3D_implicit(u,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,dt) - File Exchange

Code covered by the BSD License

Highlights from
Image Edge Enhancing Coherence Filter Toolbox

showcs3(varargin) SHOWCS3 M-file for showcs3.fig
CoherenceFilter(u,Options) This function COHERENCEFILTER will perform Anisotropic Diffusion of a
D=derivatives(I,option) Sobel like derivatives with Scharr rotation invariance stencil notations.
I=imgaussian(I,sigma,siz) IMGAUSSIAN filters an 1D, 2D color/greyscale or 3D image with an
[Dxx,Dxy,Dxz,Dyy,Dyz,Dzz]=Str... From Structure Tensor to Diffusion Tensor
[Dxx,Dxy,Dyy]=ConstructDiffus... Construct the edge preserving diffusion tensor D = [a,b;b,c] such as
[Iout]=ut_gauss(varargin) ut_gauss - 2-D filtering using Gaussian masks
[Jxx, Jxy, Jxz, Jyy, Jyz, Jzz... This function calculates the 3D Structure Tensor
[Jxx, Jxy, Jyy]=StructureTens... This function calculates the 2D
[Lambda1,Lambda2,I2x,I2y,I1x,... This function computes the eigenvectors and eigen values of the 2D image
[mu3,mu2,mu1,v3x,v3y,v3z,v2x,... This function calculates the eigen values and vectors, of the 3D image hessian.
[x,fval,exitflag,output,grad]... FMINLBFGS finds a local minimum of a function of several variables.
[x11,x12,x22]=hessian(x)
a=clmp(a,mina,maxa) Function clmp to boundaries
diffusion_scheme_3D_standard(... Standard Discretization of 3D image diffusion.
err=error_diffusion_scheme_2D... Parameters must be positive. Otherwise the kernels will go to a more
err=error_diffusion_scheme_3D... Parameters must be positive. Otherwise the kernels will go to a more
err=fitgaussian(I,a,t)
optimize_scheme_2D.m This script finds optimized kernel values for in the novel discretization
optimize_scheme_3D.m Initial values for unknowns in optimal derivative kernels
u=diffusion_scheme_2D_implici... Diffusion scheme as introduced by Weickert "Anisotropic Diffusion
u=diffusion_scheme_2D_non_neg... The Basic non_negativity diffusion equation. (Can be found in "A Scheme for
u=diffusion_scheme_2D_novel(u... par=[0.007520981141059, 0.049564810649554, 0.031509665995882, ...
u=diffusion_scheme_2D_rotatio... Most diffusion discretizations are not rotation-invariant, probably
u=diffusion_scheme_2D_standar... The standard diffusion equation. (Can be found in "A Scheme for
u=diffusion_scheme_3D_implici... This 3D Diffusion scheme is based on introduced by Weickert
u=diffusion_scheme_3D_non_neg... This is a basic non_negativity discretization of the 3D diffusion
u=diffusion_scheme_3D_novel(u...
u=diffusion_scheme_3D_novel_g...
u=diffusion_scheme_3D_rotatio... Most diffusion discretizations are not rotation-invariant. This is an
CoherenceFilterStep2D.m
CoherenceFilterStep3D.m
compile_c_files.m
View all files

from Image Edge Enhancing Coherence Filter Toolbox by Dirk-Jan Kroon
Advanced 2D/3D noise removal and edge enhancing with anisotropic diffusion filtering ( Weickert )

u=diffusion_scheme_3D_implicit(u,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,dt)

function u=diffusion_scheme_3D_implicit(u,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,dt)
% This 3D Diffusion scheme is based on introduced by Weickert  
% "Anisotropic Diffusion  in Image Processing", Thesis 1996. 
%
% !! Scheme is unstable, and not ready to use yet.
% 
% Function is written by D.Kroon University of Twente (September 2009)


%  Positive and negative image indices
[N1,N2,N3] = size(u);
id = [2:N1,N1]; iu = [1,1:N1-1];
ir = [2:N2,N2]; il = [1,1:N2-1];
ib = [2:N3,N3]; iv = [1,1:N3-1];

% Diagonal neighbor intensity updates (amount of greyvalue flow from
% neighbor) 
Dulc = max(0,-Dxy(iu,il,:)) + max(0,-Dxy); 
Ddrc = max(0,-Dxy(id,ir,:)) + max(0,-Dxy);
Ddlc = max(0, Dxy(id,il,:)) + max(0, Dxy);
Durc = max(0, Dxy(iu,ir,:)) + max(0, Dxy);

Ducv = max(0,-Dxz(iu,:,iv)) + max(0,-Dxz); 
Ddcb = max(0,-Dxz(id,:,ib)) + max(0,-Dxz);
Ddcv = max(0, Dxz(id,:,iv)) + max(0, Dxz);
Ducb = max(0, Dxz(iu,:,ib)) + max(0, Dxz);

Dclv = max(0,-Dyz(:,il,iv)) + max(0,-Dyz); 
Dcrb = max(0,-Dyz(:,ir,ib)) + max(0,-Dyz);
Dcrv = max(0, Dyz(:,ir,iv)) + max(0, Dyz);
Dclb = max(0, Dyz(:,il,ib)) + max(0, Dyz);

% Normal neighbor intensity updates 
Ducc = (Dxx(iu,:,:) + Dxx) - (abs(Dxy(iu,:,:)) + abs(Dxy)) - (abs(Dxz(iu,:,:)) + abs(Dxz));
Ddcc = (Dxx(id,:,:) + Dxx) - (abs(Dxy(id,:,:)) + abs(Dxy)) - (abs(Dxz(id,:,:)) + abs(Dxz));
Dclc = (Dyy(:,il,:) + Dyy) - (abs(Dxy(:,il,:)) + abs(Dxy)) - (abs(Dyz(:,il,:)) + abs(Dyz));
Dcrc = (Dyy(:,ir,:) + Dyy) - (abs(Dxy(:,ir,:)) + abs(Dxy)) - (abs(Dyz(:,ir,:)) + abs(Dyz));
Dccv = (Dzz(:,:,iv) + Dzz) - (abs(Dxz(:,:,iv)) + abs(Dxz)) - (abs(Dyz(:,:,iv)) + abs(Dyz));
Dccb = (Dzz(:,:,ib) + Dzz) - (abs(Dxz(:,:,ib)) + abs(Dxz)) - (abs(Dyz(:,:,ib)) + abs(Dyz));

% Normalization term to preserve region average grey value
Dccc = -(Dzz(:,:,iv) + 2*Dzz + Dzz(:,:,ib)) + ... 
       -(Dyy(:,il,:) + 2*Dyy + Dyy(:,ir,:)) + ... 
       -(Dxx(iu,:,:) + 2*Dxx + Dxx(id,:,:)) + ...  
        -max(0,-Dxy(iu,il,:)) + ...
        -max(0,-Dxy(id,ir,:)) + ...
        -max(0, Dxy(id,il,:)) + ...
        -max(0, Dxy(iu,ir,:)) + ...
        -max(0,-Dxz(iu,:,iv)) + ...
        -max(0,-Dxz(id,:,ib)) + ...
        -max(0, Dxz(id,:,iv)) + ...
        -max(0, Dxz(iu,:,ib)) + ...
        -max(0,-Dyz(:,il,iv)) + ...
        -max(0,-Dyz(:,ir,ib)) + ...
        -max(0, Dyz(:,ir,iv)) + ...
        -max(0, Dyz(:,il,ib)) + ...
       (abs(Dxy(:,il,:)) + abs(Dxy(:,ir,:)) + ...
        abs(Dxy(id,:,:)) + abs(Dxy(iu,:,:)) + ...
        abs(Dxz(id,:,:)) + abs(Dxz(iu,:,:)) + ...
        abs(Dxz(:,:,iv)) + abs(Dxz(:,:,ib)) + ...
        abs(Dyz(:,il,:)) + abs(Dyz(:,ir,:)) + ...
        abs(Dyz(:,:,iv)) + abs(Dyz(:,:,ib)) + ...
        2*abs(Dxy) + 2*abs(Dxz) + 2*abs(Dyz));

% Calculate Diffusion update 
du=0.5*dt*( Ducc.*u(iu,: ,: )  +  Ddcc.*u(id,: ,: )  + ...
            Dclc.*u(: ,il,: )  +  Dcrc.*u(: ,ir,: )  + ...
            Dccv.*u(: ,: ,iv)  +  Dccb.*u(: ,: ,ib)  + ...
            Dulc.*u(iu,il,: )  +  Ddlc.*u(id,il,: )  + ...
            Durc.*u(iu,ir,: )  +  Ddrc.*u(id,ir,: )  + ...
            Ducv.*u(iu,: ,iv)  +  Ddcb.*u(id,: ,ib)  + ...
            Ducb.*u(iu,: ,ib)  +  Ddcv.*u(id,: ,iv)  + ...
            Dclv.*u(: ,il,iv)  +  Dcrb.*u(: ,ir,ib)  + ...
            Dclb.*u(: ,il,ib)  +  Dcrv.*u(: ,ir,iv));

% Perform the edge preserving diffusion filtering on the image
u = (u + du )./ (1-dt*0.5*Dccc);  %

Highlights from Image Edge Enhancing Coherence Filter Toolbox

Highlights from
Image Edge Enhancing Coherence Filter Toolbox