Total Variation and Mean Curvature PDEs on the Homogeneous Space of Positions and Orientations

Bart (Smets^1) • Remco (Duits^1) • Etienne St - (Onge^2) • Jim (Portegies^1)

(( )^1)CASA, Eindhoven University of Technology

(( )^2)SCIL,Sherbrooke Connectivity Imaging Lab

b.m.n.smets@tue.nl
bmnsmets.com

This notebook  accompanies the paper of the same name and replicates its 2D denoising experiments (except the BM3D experiments which were done in MatLab with the code from https://www.cs.tut.fi/~foi/GCF-BM3D/). It requires the LieAnalysis package installed into Mathematica, see www.lieanalysis.nl.

Setup

$HistoryLength=1; (* Improves performance somewhat *)
<< LieAnalysis` (* Import LieAnalysis package *)

Images used in the paper.

collagenImage = ;
monaLisaImage = ;
spiralImage = ;

Select grayscale image.

sourceImage = ColorConvert[monaLisaImage, "Grayscale"];

Generate noise.

noise = RandomVariate[NormalDistribution[0,0.2], Dimensions[ImageData[sourceImage]]];
noise = GaussianFilter[noise, 2]; (* Comment this line out to turn off noise correlation *)

noisyImage = Image[Clip[ImageData[sourceImage]+noise], ImageSize->Medium]

Metric constants.

Ds=1;
Da=1/100;
eGothic = 1;

Orientation Score

Generate orientation score (with 8 effective orientations).

(*Orientations*)
Nθ=8; 
(*Angular resolution*)
sθ = π/Nθ;

os = OrientationScoreTransform[noisyImage, "Orientations"-> 2*Nθ, "WaveletSize"-> 75];

Local Adaptive Frames

Helper functions.

TensorToGaugeFrame=Compile[{{tensor,_Real,2},{gaugeFrame,_Real,2}},gaugeFrame\[Transpose].tensor.gaugeFrame,RuntimeAttributes->{Listable}];

CalcGaugeFrames [os_, Ds_, Da_, κ_, sθ_, σs_:0, σa_:0] := Module[{A,H, M,MHM,B1, Chi, Nu, B2, B3},
H = OrientationScoreTensor[os, {1, sθ}, "Hessian"]["Data"];
H = H/Max[Abs[H]];
If[{σs,σa} =!= {0,0},
H = LieAnalysis`Common`ExternalRegularization[H, os["FullOrientationList"], {σs,σa}];
];
If[κ==0
,
H = H[[All,All,{1,3},{1,3}]];
M = DiagonalMatrix[{Sqrt[Ds], Sqrt[Da]}];
,
M = DiagonalMatrix[{Sqrt[Ds],κ Sqrt[Ds],  Sqrt[Da]}];
];
(*MHM = Map[M.#.M & ,H, {3}];
A = Map[#\[Transpose].#&,MHM,{3}];*)
A = Map[(M^2).#\[Transpose].(M^2).# &, H, {3}];

B1 = Map[(Eigenvectors[#, Method->"Direct"][[3]]) &, A, {3}];
B1 = Map[Sign[#[[1]]] #/Sqrt[#[[1]]^2/Ds+ #[[2]]^2/(κ^2 Ds)+#[[3]]^2/Da]&,B1,{3}];

Chi = Mod[ArcTan[B1[[All,All,All,1]], B1[[All,All,All,2]]],Pi,-Pi/2]; Nu = ArcTan[Sqrt[B1[[All,All,All,1]]^2 + B1[[All,All,All,2]]^2], B1[[All,All,All,3]] Sign[B1[[All,All,All,1]]]];

B2 = MapThread[ {-Sin[#1], Cos[#1], 0} &,{Chi, Nu}, 3];
B2 = Map[#/Sqrt[#[[1]]^2/Ds+ #[[2]]^2/(κ^2 Ds)+#[[3]]^2/Da]&,B2,{3}];

B3 = MapThread[{Cos[#1]*Sin[#2], Sin[#2]*Sin[#1], Cos[#2]} & , {Chi, Nu}, 3];
B3 = Map[#/Sqrt[#[[1]]^2/Ds+ #[[2]]^2/(κ^2 Ds)+#[[3]]^2/Da]&,B3,{3}];

Return[Transpose[{B1,B2,B3}, {4,1,2,3,5}]]
];

Fit gauge frames.

gaugeFrames = CalcGaugeFrames[os, Ds, Da, 1, sθ, 6., 4.*sθ];

Isotropy Factor

hessian = OrientationScoreTensor[os, {3./2, 3./2*sθ}, "Hessian"];
hessianGauge = TensorToGaugeFrame[hessian["Data"], gaugeFrames];

c=0.2;

orConfG = Abs@(hessianGauge[[All,All,All,2,2]]+hessianGauge[[All,All,All,3,3]]);
orConfL = Abs@(hessian["Data"][[All,All,All,2,2]]+hessian["Data"][[All,All,All,3,3]]);
αG = Exp[-Rescale[orConfG]/c];
αL = Exp[-Rescale[orConfL]/c];

PDEs

Iterative denoisers which returns the denoised image and L1 and L2 error series.

Left Invariant MCF

DoLeftMCF[os_,α_, t_, dt_, Ds_,Da_, source_, κ_:1] := Module[
{out = os, progress=0.0, diff, grad, grad1, grad2, grad3, gradNorm, gradGN1, gradGN2, gradGN3, gradGN, L1, L2, img = InverseOrientationScoreTransform[os]},
Monitor[
L1 = {ImageDistance[source, img, DistanceFunction->ManhattanDistance]};
L2 = {ImageDistance[source, img, DistanceFunction->EuclideanDistance]};
Do[
grad = OrientationScoreTensor[out, "Gradient"]["Data"];
grad1 = Ds*grad[[All,All,All,1]];
grad2 = Ds*κ^2*grad[[All,All,All,2]];
grad3 = Da*grad[[All,All,All,3]];
gradNorm = Sqrt[grad1^2 / Ds+grad2^2/(Ds*κ^2)+grad3^2/Da+0.001^2];

gradNorm = GaussianFilter[gradNorm, {{1,1,0}}, Padding->"Fixed"];
gradNorm = GaussianFilter[gradNorm, {{0,0,1/2}}, Padding-> "Periodic"];

gradGN1 = OrientationScoreTensor[Affix[out, "Data"-> grad1/gradNorm], "Gradient"]["Data"];
gradGN2 = OrientationScoreTensor[Affix[out, "Data"-> grad2/gradNorm], "Gradient"]["Data"];
gradGN3 = OrientationScoreTensor[Affix[out, "Data"-> grad3/gradNorm], "Gradient"]["Data"];

gradGN = MapThread[{#1,#2,#3}&,{gradGN1, gradGN2, gradGN3},4];

diff = gradNorm*(gradGN[[All,All,All,1,1]]
+α*gradGN[[All,All,All,2,2]]
+α*gradGN[[All,All,All,3,3]]);

out = Affix[out, "Data"-> out["Data"] + dt*diff];
img = InverseOrientationScoreTransform[out];
AppendTo[L1, ImageDistance[source, img, DistanceFunction->ManhattanDistance]];
AppendTo[L2, ImageDistance[source, img, DistanceFunction->EuclideanDistance]];
progress += dt/t;
, Round[t/dt]
];
, {ProgressIndicator[progress], 
Image[img, ImageSize->800]}];
Return[{out, L1, L2}];
];

Left Invariant TVF

DoLeftTVF[os_,α_, t_, dt_, Ds_, Da_, source_, κ_:1] := Module[
{out = os, progress=0.0, diff, grad, grad1, grad2, grad3, gradNorm, gradGN1, gradGN2, gradGN3, gradGN,λ=50, L1, L2, img=InverseOrientationScoreTransform[os]},
Monitor[
L1 = {ImageDistance[source, img, DistanceFunction->ManhattanDistance]};
L2 = {ImageDistance[source, img, DistanceFunction->EuclideanDistance]};
out = Affix[out, "Data"-> λ * out["Data"]];
Do[
grad = OrientationScoreTensor[out, "Gradient"]["Data"];
grad1 = Ds*grad[[All,All,All,1]];
grad2 = Ds*κ^2*grad[[All,All,All,2]];
grad3 = Da*grad[[All,All,All,3]];
gradNorm = Sqrt[grad1^2/Ds+grad2^2/(Ds*κ^2)+grad3^2/Da+0.001^2];

gradNorm = GaussianFilter[gradNorm, {{1,1,0}}, Padding->"Fixed"];
gradNorm = GaussianFilter[gradNorm, {{0,0,1/2}}, Padding-> "Periodic"];

gradGN1 = OrientationScoreTensor[Affix[out, "Data"-> grad1/gradNorm], "Gradient"]["Data"];
gradGN2 = OrientationScoreTensor[Affix[out, "Data"-> grad2/gradNorm], "Gradient"]["Data"];
gradGN3 = OrientationScoreTensor[Affix[out, "Data"-> grad3/gradNorm], "Gradient"]["Data"];

gradGN = MapThread[{#1,#2,#3}&,{gradGN1, gradGN2, gradGN3},4];

diff = (gradGN[[All,All,All,1,1]]
+α*gradGN[[All,All,All,2,2]]
+α*gradGN[[All,All,All,3,3]]);

out = Affix[out, "Data"-> out["Data"] + dt*diff];
img = InverseOrientationScoreTransform[Affix[out, "Data"-> (out["Data"]/λ)]];
AppendTo[L1, ImageDistance[source, img, DistanceFunction->ManhattanDistance]];
AppendTo[L2, ImageDistance[source, img, DistanceFunction->EuclideanDistance]];
progress += dt/t;
, Round[t/dt]
];
, {ProgressIndicator[progress],
Image[img, ImageSize->800]}];
Return[{Affix[out, "Data"-> (out["Data"]/λ)],L1,L2}];
];

Gauge Frame MCF

DoGaugeMCF[os_, gaugeFrames_, α_, t_, dt_, Ds_, Da_, source_, κ_:1] := Module[
{out = os, progress=0.0, diff, grad, grad1, grad2, grad3, gradNorm, gradGN1, gradGN2, gradGN3, gradGN, gradGNGauge, L1, L2, img = InverseOrientationScoreTransform[os]},
Monitor[
L1 = {ImageDistance[source, img, DistanceFunction->ManhattanDistance]};
L2 = {ImageDistance[source, img, DistanceFunction->EuclideanDistance]};
Do[
grad = OrientationScoreTensor[out, "Gradient"]["Data"];
grad1 = Ds*grad[[All,All,All,1]];
grad2 = Ds*κ^2*grad[[All,All,All,2]];
grad3 = Da*grad[[All,All,All,3]];
gradNorm = Sqrt[grad1^2/Ds+grad2^2/(Ds*κ^2)+grad3^2/Da+0.001^2];

gradNorm = GaussianFilter[gradNorm, {{1,1,0}}, Padding->"Fixed"];
gradNorm = GaussianFilter[gradNorm, {{0,0,1/2}}, Padding-> "Periodic"];

gradGN1 = OrientationScoreTensor[Affix[out, "Data"-> grad1/gradNorm], "Gradient"]["Data"];
gradGN2 = OrientationScoreTensor[Affix[out, "Data"-> grad2/gradNorm], "Gradient"]["Data"];
gradGN3 = OrientationScoreTensor[Affix[out, "Data"-> grad3/gradNorm], "Gradient"]["Data"];

gradGN = MapThread[{#1,#2,#3}&,{gradGN1, gradGN2, gradGN3},4];
gradGNGauge = TensorToGaugeFrame[gradGN, gaugeFrames];

diff = gradNorm*(gradGNGauge[[All,All,All,1,1]]
+α*gradGNGauge[[All,All,All,2,2]]
+α*gradGNGauge[[All,All,All,3,3]]);

out = Affix[out, "Data"-> out["Data"] + dt*diff];
img = InverseOrientationScoreTransform[out];
AppendTo[L1, ImageDistance[source, img, DistanceFunction->ManhattanDistance]];
AppendTo[L2, ImageDistance[source, img, DistanceFunction->EuclideanDistance]];
progress += dt/t;
, Round[t/dt]
];
, {ProgressIndicator[progress], 
Image[img, ImageSize->800]}];
Return[{out,L1,L2}];
];

Gauge Frame TVF

DoGaugeTVF[os_, gaugeFrames_, α_, t_, dt_, Ds_, Da_, source_, κ_:1] := Module[
{out = os, progress=0.0, diff, grad, grad1, grad2, grad3, gradNorm, gradGN1, gradGN2, gradGN3, gradGN, gradGNGauge, λ=50, L1, L2, img=InverseOrientationScoreTransform[os]},
Monitor[
L1 = {ImageDistance[source, img, DistanceFunction->ManhattanDistance]};
L2 = {ImageDistance[source, img, DistanceFunction->EuclideanDistance]};
out = Affix[out, "Data"-> λ * out["Data"]];
Do[
grad = OrientationScoreTensor[out, "Gradient"]["Data"];
grad1 = Ds*grad[[All,All,All,1]];
grad2 = Ds*κ^2*grad[[All,All,All,2]];
grad3 = Da*grad[[All,All,All,3]];
gradNorm = Sqrt[grad1^2/Ds+grad2^2/(Ds*κ^2)+grad3^2/Da+0.001^2];

gradNorm = GaussianFilter[gradNorm, {{1,1,0}}, Padding->"Fixed"];
gradNorm = GaussianFilter[gradNorm, {{0,0,1/2}}, Padding-> "Periodic"];
gradGN1 = OrientationScoreTensor[Affix[out, "Data"-> grad1/gradNorm], "Gradient"]["Data"];
gradGN2 = OrientationScoreTensor[Affix[out, "Data"-> grad2/gradNorm], "Gradient"]["Data"];
gradGN3 = OrientationScoreTensor[Affix[out, "Data"-> grad3/gradNorm], "Gradient"]["Data"];

gradGN = MapThread[{#1,#2,#3}&,{gradGN1, gradGN2, gradGN3},4];
gradGNGauge = TensorToGaugeFrame[gradGN, gaugeFrames];

diff = (gradGNGauge[[All,All,All,1,1]]
+α*gradGNGauge[[All,All,All,2,2]]
+α*gradGNGauge[[All,All,All,3,3]]);

out = Affix[out, "Data"-> out["Data"] + dt*diff];
img = InverseOrientationScoreTransform[Affix[out, "Data"-> (out["Data"]/λ)]];
AppendTo[L1, ImageDistance[source, img, DistanceFunction->ManhattanDistance]];
AppendTo[L2, ImageDistance[source, img, DistanceFunction->EuclideanDistance]];
progress += dt/t;
, Round[t/dt]
];
, {ProgressIndicator[progress],
Image[img, ImageSize->800]}];
Return[{Affix[out, "Data"-> (out["Data"]/λ)], L1, L2}];
];

Perona - Malik

DoPM[os_,t_,dt_,source_] := Module[
{out = InverseOrientationScoreTransform[os], progress=0.0, L1, L2, img=InverseOrientationScoreTransform[os]},
Monitor[
L1 = {ImageDistance[source, img, DistanceFunction->ManhattanDistance]};
L2 = {ImageDistance[source, img, DistanceFunction->EuclideanDistance]};
Do[

out = PeronaMalikFilter[out,dt];
AppendTo[L1, ImageDistance[source, out, DistanceFunction->ManhattanDistance]];
AppendTo[L2, ImageDistance[source, out, DistanceFunction->EuclideanDistance]];
progress += dt/t;
, Round[t/dt]
];
, {ProgressIndicator[progress],
Image[out, ImageSize->800]}];
Return[{out, L1, L2}];
];

Denoising

(* Turn coherence enhancement on/off *)
coherenceEnhancement = True;

t=6.; (* Processing time *)
dt = 0.1; (* processing step *)

ts = Range[0,t,dt];
nullImage = 0*sourceImage;
sourceL1 = ImageDistance[sourceImage, nullImage, DistanceFunction-> ManhattanDistance];
sourceL2 = ImageDistance[sourceImage, nullImage, DistanceFunction-> EuclideanDistance];

If[coherenceEnhancement, αinG = αG; αinL = αL, αinG = 1; αinL = 1];

Run experiments.

{pmImage, pm1, pm2} = DoPM[os, 10 t, 10 dt, sourceImage]; (* bring PM timescale into line *)
{lmImage, lm1, lm2} = DoLeftMCF[os,αinL, t, dt, Ds, Da, sourceImage, eGothic];
{gmImage, gm1, gm2} = DoGaugeMCF[os, gaugeFrames,αinG, t, dt, Ds, Da, sourceImage, eGothic];
{ltImage, lt1, lt2} = DoLeftTVF[os,αinL, t, dt, Ds, Da, sourceImage, eGothic];
{gtImage, gt1, gt2} = DoGaugeTVF[os, gaugeFrames, αinG, t, dt, Ds, Da, sourceImage, eGothic];

Comparison

L1/L2 Errors

colorPM = Red;
colorMCF = Lighter@Blue;
colorTVF = Darker@Green;

pm1plot = Style[Transpose[{ts, pm1/sourceL1}], colorPM];

lm1plot = Style[Transpose[{ts, lm1/sourceL1}], colorMCF, Dashed];
lt1plot = Style[Transpose[{ts, lt1/sourceL1}],  colorTVF,Dashed];

gm1plot = Style[Transpose[{ts, gm1/sourceL1}], colorMCF];
gt1plot = Style[Transpose[{ts, gt1/sourceL1}], colorTVF];

pm2plot = Style[Transpose[{ts, pm2/sourceL2}], colorPM];

lm2plot = Style[Transpose[{ts, lm2/sourceL2}],  colorMCF,Dashed];
lt2plot = Style[Transpose[{ts, lt2/sourceL2}], colorTVF,Dashed];

gm2plot = Style[Transpose[{ts, gm2/sourceL2}],colorMCF ];
gt2plot = Style[Transpose[{ts, gt2/sourceL2}], colorTVF];

leg=LineLegend[{colorPM, Directive[colorMCF, Dashed],colorMCF,Directive[colorTVF, Dashed],  colorTVF},{"Perona-Malik 2D", "Left Invariant MCF","Gauge MCF", "Left Invariant TVF",  "Gauge TVF"}, LegendLayout->{"Row",1}];

Legended[Row[{
ListLinePlot[{pm1plot, lm1plot, lt1plot, gm1plot, gt1plot}, 
PlotLabel->Text[Style[TraditionalForm["(Subscript[L, 1]) Errors"],Large]],
ImageSize->440
]
,
ListLinePlot[{pm2plot, lm2plot, lt2plot, gm2plot, gt2plot}, 
PlotLabel->Text[Style[TraditionalForm["(Subscript[L, 2]) Errors"],Large]],
ImageSize->440
]
},Spacer[40], ImageSize->Full],  Placed[leg, Below]]

Time to minimal L2 error

tPM = dt*(FirstPosition[pm2, Min[pm2]][[1]]-1)
tLM = dt*(FirstPosition[lm2, Min[lm2]][[1]]-1)
tGM = dt*(FirstPosition[gm2, Min[gm2]][[1]]-1)
tLT = dt*(FirstPosition[lt2, Min[lt2]][[1]]-1)
tGT = dt*(FirstPosition[gt2, Min[gt2]][[1]]-1)

Qualitative results

GraphicsGrid[{{sourceImage, noisyImage},
{pmImage, InverseOrientationScoreTransform@os},
{InverseOrientationScoreTransform@lmImage, InverseOrientationScoreTransform@gmImage},
{InverseOrientationScoreTransform@ltImage, InverseOrientationScoreTransform@gtImage}}, ImageSize->Full, ItemAspectRatio->ImageAspectRatio[sourceImage]]

Oversmooth an image a set amount past the optimal L2 point.

toomuch=1.5;
{qgtImage, qgt1, qgt2} = DoGaugeTVF[os, gaugeFrames, αinG, toomuch*tGT, dt, Ds, Da, sourceImage, eGothic];

InverseOrientationScoreTransform@qgtImage

PSNR

PSNR[img_, source_] := Module[{w=ImageDimensions[img][[1]], h=ImageDimensions[img][[2]]},
-10*Log10[Total[(ImageData[img]- ImageData[source])^2,2]/(h w)]
];

PSNR[noisyImage, sourceImage]
PSNR[pmImage, sourceImage]
PSNR[InverseOrientationScoreTransform[lmImage],sourceImage]
PSNR[InverseOrientationScoreTransform[gmImage],sourceImage]
PSNR[InverseOrientationScoreTransform[ltImage],sourceImage]
PSNR[InverseOrientationScoreTransform[gtImage],sourceImage]