Image Edge Enhancing Coherence Filter Toolbox: optimize_scheme_3D.m - 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 )

optimize_scheme_3D.m

functionname='optimize_scheme_3D.m';
functiondir=which(functionname);
functiondir=functiondir(1:end-length(functionname));
addpath([functiondir '../functions2D'])
addpath([functiondir '../functions3D'])
addpath([functiondir '../functions'])

%% Initial values for unknowns in optimal derivative kernels
par =[0.029180542350452,0.051009964656270,0.056702441663829,0.000006816925755, ...
          0.000008999487006,0.000009650996104,0.004751977565242,0.000082895902805, ...
          0.000000310453682,0.002474067467941,0.000007941958352,0.004974632072065, ...
          0.000174128474993,0.004854111687777,0.090614087403327,0.000116874531236, ...
          0.005034186669466,0.000373963590204,0.017456316783794,0.000008810909572, ...
          0.194390591201552,0.000005600928728,0.000003816386672,0.000003302536833, ...
          0.000014736927086,0.000000132136404,0.000005381620877,0.000003078916993, ...
          0.000005671182049,0.000004345622765,0.000004236931951,0.000002758996948, ...
          0.043827391390423];
      
%% Make a test image with, circles of different spatial frequencies
n=60; a=8;  
[x,y,z]=ndgrid(-n:n,-n:n,-n:n);
x=1.1*x/(n/a); 
y=1.1*y/(n/a);  
z=1.1*z/(n/a);
r=sqrt(sqrt(x.^2+y.^2+z.^2))*(n/a)^2;
u_perfect=sin(r);

u_noise=u_perfect+rand(size(u_perfect))*0.05-0.025;
u_perfect=single(u_perfect);
u_noise=single(u_noise);
%% Optimization options
options_opt.dt=0.5;
options_opt.alpha=50000;
options_opt.itt1=30;
options_opt.itt2=5;
options_opt.verbose=false;

%% Make the Diffusion kernel eigenvectors
Options=struct('sigma', 1, 'rho', 1, 'TensorType', 1, 'eigenmode',0,'C', 1e-10, 'm',1,'alpha',0.001,'lambda_e',0.02,'lambda_c',0.02,'lambda_h',0.5,'RealDerivatives',false);

% Gaussian smooth the image, for better gradients
usigma=imgaussian(u_perfect,Options.sigma,4*Options.sigma);

% Calculate the gradients
ux=derivatives(usigma,'x');
uy=derivatives(usigma,'y');
uz=derivatives(usigma,'z');

% Compute the 3D structure tensors J of the image
[Jxx, Jxy, Jxz, Jyy, Jyz, Jzz] = StructureTensor3D(ux,uy,uz, Options.rho);

% Gradient magnitude squared
gradA=ux.^2+uy.^2+uz.^2;

% Free memory
clear ux; clear uy; clear uz;

% Compute the eigenvectors and eigenvalues of the hessian and directly
% use the equation of Weickert to convert them to diffusion tensors
[Dxx,Dxy,Dxz,Dyy,Dyz,Dzz]=StructureTensor2DiffusionTensor3D(Jxx,Jxy,Jxz,Jyy,Jyz,Jzz,gradA,Options); 

% Free memory
clear J*;

%% Find the optimal derivative values
for i=1:10,
    par= fminlbfgs (@(x)error_diffusion_scheme_3D(u_perfect,u_noise,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,x,options_opt),par,struct('Display','iter','MaxIter',1000,'TolX',1e-6));
    par= fminsearch(@(x)error_diffusion_scheme_3D(u_perfect,u_noise,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,x,options_opt),par,struct('Display','iter','MaxIter',1000,'TolX',1e-6));
    save('par3d');
end

%% Show the results
error_diffusion_scheme_3D(u_perfect,u_noise,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,par,options_opt)

Highlights from Image Edge Enhancing Coherence Filter Toolbox

Highlights from
Image Edge Enhancing Coherence Filter Toolbox