energy function calculation for TV models

4 files changed, 81 insertions, 2 deletions

diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
index 973d060..84889d7 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -21,6 +21,7 @@ fprintf('Denoise a volume using the ROF-TV model (CPU) \n');
 tau_rof = 0.0025; % time-marching constant 
 iter_rof = 300; % number of ROF iterations
 tic; u_rof = ROF_TV(single(vol3D), lambda_reg, iter_rof, tau_rof); toc; 
+energyfunc_val_rof = TV_energy(single(u_rof),single(vol3D),lambda_reg);  % get energy function value
 figure; imshow(u_rof(:,:,15), [0 1]); title('ROF-TV denoised volume (CPU)');
 %%
 % fprintf('Denoise a volume using the ROF-TV model (GPU) \n');
@@ -33,6 +34,7 @@ fprintf('Denoise a volume using the FGP-TV model (CPU) \n');
 iter_fgp = 300; % number of FGP iterations
 epsil_tol =  1.0e-05; % tolerance
 tic; u_fgp = FGP_TV(single(vol3D), lambda_reg, iter_fgp, epsil_tol); toc; 
+energyfunc_val_fgp = TV_energy(single(u_fgp),single(vol3D),lambda_reg); % get energy function value
 figure; imshow(u_fgp(:,:,15), [0 1]); title('FGP-TV denoised volume (CPU)');
 %%
 % fprintf('Denoise a volume using the FGP-TV model (GPU) \n');
@@ -45,6 +47,7 @@ fprintf('Denoise a volume using the SB-TV model (CPU) \n');
 iter_sb = 150; % number of SB iterations
 epsil_tol =  1.0e-05; % tolerance
 tic; u_sb = SB_TV(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; 
+energyfunc_val_sb = TV_energy(single(u_sb),single(vol3D),lambda_reg);  % get energy function value
 figure; imshow(u_sb(:,:,15), [0 1]); title('SB-TV denoised volume (CPU)');
 %%
 % fprintf('Denoise a volume using the SB-TV model (GPU) \n');
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 4a0a19a..526d21c 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -12,13 +12,14 @@ lambda_reg = 0.03; % regularsation parameter for all methods
 %%
 fprintf('Denoise using the ROF-TV model (CPU) \n');
 tau_rof = 0.0025; % time-marching constant 
-iter_rof = 2000; % number of ROF iterations
+iter_rof = 750; % number of ROF iterations
 tic; u_rof = ROF_TV(single(u0), lambda_reg, iter_rof, tau_rof); toc; 
+energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg);  % get energy function value
 figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)');
 %%
 % fprintf('Denoise using the ROF-TV model (GPU) \n');
 % tau_rof = 0.0025; % time-marching constant 
-% iter_rof = 2000; % number of ROF iterations
+% iter_rof = 750; % number of ROF iterations
 % tic; u_rofG = ROF_TV_GPU(single(u0), lambda_reg, iter_rof, tau_rof); toc;
 % figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)');
 %%
@@ -26,7 +27,9 @@ fprintf('Denoise using the FGP-TV model (CPU) \n');
 iter_fgp = 1000; % number of FGP iterations
 epsil_tol =  1.0e-06; % tolerance
 tic; u_fgp = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; 
+energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg); % get energy function value
 figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)');
+
 %%
 % fprintf('Denoise using the FGP-TV model (GPU) \n');
 % iter_fgp = 1000; % number of FGP iterations
@@ -38,6 +41,7 @@ fprintf('Denoise using the SB-TV model (CPU) \n');
 iter_sb = 150; % number of SB iterations
 epsil_tol =  1.0e-06; % tolerance
 tic; u_sb = SB_TV(single(u0), lambda_reg, iter_sb, epsil_tol); toc; 
+energyfunc_val_sb = TV_energy(single(u_sb),single(u0),lambda_reg);  % get energy function value
 figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)');
 %%
 % fprintf('Denoise using the SB-TV model (GPU) \n');
diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
index ec799bd..a445e99 100644
--- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
@@ -23,6 +23,9 @@ movefile TNV.mex* ../installed/
 mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 movefile NonlDiff.mex* ../installed/
 
+mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile TV_energy.mex* ../installed/
+
 delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h
 
 fprintf('%s \n', 'All successfully compiled!');
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c
new file mode 100644
index 0000000..421bd4c
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c
@@ -0,0 +1,69 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "utils.h"
+/*
+ * Function to calculate TV energy value with respect to the denoising variational problem
+ * 
+ * Input:
+ * 1. Denoised Image/volume
+ * 2. Original (noisy) Image/volume
+ * 3. lambda - regularisation parameter 
+ * 
+ * Output:
+ * 1. Energy function value
+ * 
+ */
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, dimX, dimY, dimZ;
+    const int  *dim_array;
+    float *Input, *Input0, lambda;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs != 3)) mexErrMsgTxt("3 inputs: Two images or volumes of the same size required, estimated and the original (noisy), regularisation parameter");
+    
+    Input  = (float *) mxGetData(prhs[0]); /* Denoised Image/volume */
+    Input0  = (float *) mxGetData(prhs[1]); /* Original (noisy) Image/volume */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularisation parameter */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    /*output energy function value */
+    plhs[0] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
+    float *funcvalA = (float *) mxGetData(plhs[0]);
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+		TV_energy2D(Input, Input0, funcvalA, lambda, dimX, dimY);
+		}
+    if (number_of_dims == 3) {
+        TV_energy3D(Input, Input0, funcvalA, lambda, dimX, dimY, dimZ);
+    }
+}