The test image

image=ImageResize[ColorConvert[ExampleData[{"TestImage","Elaine"}],"Grayscale"],256];
f0=ImageData[image];
f=ImageData[ImageAdjust[ImageEffect[Image[f0],{"GaussianNoise",0.1}]]];
GraphicsRow[{Image[f0],Image[f]},ImageSize->800]

Linear isotropic diffusion

These function compute the second order derivatives as finite differences along the array axes

Unprotect[xx,yy];xx=.;yy=.;Protect[xx,yy];
(*Second order central differences*)
(Subscript[f_, xx]):=RotateLeft[f,{1,0}]-2f+RotateLeft[f,{-1,0}]
(Subscript[f_, yy]):=RotateLeft[f,{0,1}]-2f+RotateLeft[f,{0,-1}]

Iteratively solve the diffusion PDE

(*Evolution settings*)
dt=0.01;(*Step size in time; Lower this for more a slower evolution*)
tFinal=10;(*Final time of the evolution*)
(*Initialize solution*)
u=u0=f;(*Initialize u at t=0*)
(*This print statement dynamically displays the current solution*)
PrintTemporary[Dynamic[Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[t],ImageSize->400]]];
(*Do the time evolution*)
Do[
(*The update step*)
u=u+dt((Subscript[u, xx])+(Subscript[u, yy]))
,{t,0,tFinal,dt}];
(*Plot the final result*)
Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[tFinal],ImageSize->400]

Non-linear isotropic diffusion (Perona-Malik diffusion)

Introduce a diffusion constant based on the gradient (low diffusion at edges)

(*Parameters*)
σ=1;(*Scale/standard dev. at which the (Gaussian) gradient magnitude is computed*)
(*Compute the gradient magnitude*)
gradMagn=GradientFilter[f,{3σ,σ}];
K=Quantile[Flatten[gradMagn],0.5];(*The Perona-Malik constant (lower means the less diffusion at edges). A good choice is the 50% quantile of the gradient magnitude map*)
c=(1/(1+((((gradMagn/K)))^2)));
(*Show the initial settings*)
GraphicsRow[{Show[Image[f],PlotLabel->"Input Image"],Show[ImageAdjust[Image[gradMagn]],PlotLabel->"Gradient magnitude"],Show[Image[c],PlotLabel->"Diffusion constant"]},ImageSize->800]

Iteratively solve the diffusion PDE

(*Evolution settings*)
dt=0.01;(*Step size in time; Lower this for more a slower evolution*)
tFinal=10;(*Final time of the evolution*)
(*Initialize solution*)
u=u0=f;(*Initialize u at t=0*)
(*This print statement dynamically displays the current solution*)
PrintTemporary[Dynamic[Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[t],ImageSize->400]]];
(*Do the time evolution*)
Do[
(*The update step, now with diffusion constant c*)
u=u+dt c((Subscript[u, xx])+(Subscript[u, yy]))
,{t,0,tFinal,dt}];
(*Plot the final result*)
Show[ImageAdjust[Image[u]],PlotLabel->"t="<>ToString[tFinal],ImageSize->400]

Anisotropic linear diffusion (fixed diffusion matrix for each (x,y))

Now we write diffusion in divergence form, and now we need a function to compute first order derivatives

Unprotect[xf,yf,xb,yb];Remove[xf,yf,xb,yb];Protect[xf,yf,xb,yb];
(*Forward differences*)
(Subscript[f_, xf]):=RotateLeft[f,{1,0}]-f
(Subscript[f_, yf]):=RotateLeft[f,{0,1}]-f
(*Backward differences*)
(Subscript[f_, xb]):=f-RotateLeft[f,{-1,0}]
(Subscript[f_, yb]):=f-RotateLeft[f,{0,-1}]
grad[f_]:={(Subscript[f, xf]),(Subscript[f, yf])}
div[ff_]:=(Subscript[ff[[1]], xb])+(Subscript[ff[[2]], yb])

Introduce diffusion matrix

θ=π/4.;(*Apply more diffusion in the direction of theta*)
(*The diffusion matrix will have the following eigenvectors*)
e1=(({{Cos[θ]}, {Sin[θ]}}));e2=(({{-Sin[θ]}, {Cos[θ]}}));
(*And the following eigen values*)
λ1=1;
λ2=0.1;
(*The diffusion matrix*)
Dmat=λ1*e1.e1\[Transpose]+λ2*e2.e2\[Transpose];
MatrixForm[Dmat]

Solve PDE (on delta spike and the image)

(*Evolution settings*)
dt=0.025;(*Step size in time; Lower this for more a slower evolution*)
tFinal=10;(*Final time of the evolution*)
(*Initialize solution*)
u=u0=f;(*Initialize u at t=0*)
v=v0=CenterArray[1,Dimensions[f]];(*Now we also keep track of the evolution on a delta spike*)
(*This print statement dynamically displays the current solution*)
PrintTemporary[Dynamic[Show[GraphicsRow[{ImageAdjust@Image[u],ImageAdjust@Image[v]},ImageSize->800],PlotLabel->"t="<>ToString[t]]]];
(*Do the time evolution*)
Do[
(*The update step, now in divergence form with diffusion matrix "Dmat"*)
u=u+dt div[Dmat.grad[u]];
v=v+dt div[Dmat.grad[v]];
,{t,0,tFinal,dt}];
(*Plot the final result*)
Show[GraphicsRow[{ImageAdjust@Image[u],ImageAdjust@Image[v]},ImageSize->800],PlotLabel->"t="<>ToString[tFinal]]

Non-linear anisotropic diffusion (coherence enhancing diffusion)

A diffusion matrix per pixel

(Subscript[grad, σ_])[f_]:=Transpose[grad[GaussianFilter[f,{3σ,σ}]],{3,1,2}](*a gradient vector at each pixel*)
(Subscript[S, σ_,ρ_])[f_]:=GaussianFilter[Map[Outer[Times,#,#]&,(Subscript[grad, σ])[f],{2}],{{3ρ,3ρ,1,1},{ρ,ρ,0,0}}](*the structure tensor at each pixel (each component blurred spatially)*)
CoherenceFromEigenvalues[eigvals_]:=Sum[Sum[(((eigvals[[i]]-eigvals[[j]]))^2),{j,i+1,Length[eigvals]}],{i,1,Length[eigvals]-1}]
StructTensToDmat[S_,constC_,α_]:=With[{eigsys=N@Eigensystem[S]},With[{κ=CoherenceFromEigenvalues[eigsys[[1]]]},With[{λm=If[κ==0.,α,α+(1-α)Exp[((-constC)/κ)]]},
Sum[Outer[Times,eigsys[[2,i]],eigsys[[2,i]]]*α,{i,1,Length[eigsys[[1]]]-1}]+Outer[Times,eigsys[[2,-1]],eigsys[[2,-1]]]*λm
]]]

gradvec[f_]:=Transpose[{grad[f]},{4,3,1,2}](*a column vector at each pixel*)
divvec[ff_]:=(Subscript[ff[[All,All,1,1]], xb])+(Subscript[ff[[All,All,2,1]], yb])

Construct diffusion matrix field

constC=0.0000001;
constC=0.000000001;(*A high value means only blurring in regions which are really anisotropic, low value means blurring every where*)
α=0.05;(*Anisotropy of the blurring along the preferred direction (the width of van Gogh-style pencil stroke)*)
σ=1.;
ρ=4.;
Dfield=Map[StructTensToDmat[#,constC,α]&,(Subscript[S, σ,ρ])[f],{2}];//AbsoluteTiming

Dfield[[10,10]]//MatrixForm

Solve PDE

(*Evolution settings*)
dt=0.05;(*Step size in time; Lower this for more a slower evolution*)
tFinal=10;(*Final time of the evolution*)
(*Initialize solution*)
u=u0=f;(*Initialize u at t=0*)
(*This print statement dynamically displays the current solution*)
PrintTemporary[Dynamic[Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[t],ImageSize->400]]];
(*Do the time evolution*)
Do[
(*The update step, now in divergence form with diffusion matrix "Dmat"*)
u=u+dt divvec[MapThread[#1.#2&,{Dfield,gradvec[u]},2]];
,{t,0,tFinal,dt}];
(*Plot the final result*)
Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[tFinal],ImageSize->400]

TV denoising

Solve PDE (on delta spike and the image)

(*Evolution settings*)
dt=0.0001;(*Step size in time; Lower this for more a slower evolution*)
tFinal=0.1;(*Final time of the evolution*)
ϵ=0.000000000000001;
(*Initialize solution*)
u=u0=f;(*Initialize u at t=0*)
(*This print statement dynamically displays the current solution*)
PrintTemporary[Dynamic[Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[t],ImageSize->400]]];
(*Do the time evolution*)
Do[
(*The update step, now in divergence form with diffusion matrix "Dmat"*)
u=u+dt div[{((#[[1]])/(ϵ+(Sqrt[((#[[1]])^2)+((#[[2]])^2)]))),((#[[2]])/(ϵ+(Sqrt[((#[[1]])^2)+((#[[2]])^2)])))}&@grad[u]];
,{t,0,tFinal,dt}];
(*Plot the final result*)
Show[ImageAdjust@Image[u],PlotLabel->"t="<>ToString[tFinal],ImageSize->400]