energy function calculation for TV models

6 files changed, 111 insertions, 2 deletions

diff --git a/Core/regularisers_CPU/utils.c b/Core/regularisers_CPU/utils.c
index a141cf4..f21d383 100644
--- a/Core/regularisers_CPU/utils.c
+++ b/Core/regularisers_CPU/utils.c
@@ -62,3 +62,32 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int
 		E_val[0] = E_Grad + E_Data;
 		return *E_val;
 }
+
+float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY, int dimZ)
+{
+	int i, j, k, i1, j1, k1, index;
+	float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f;
+	
+	/* first calculate \grad U_xy*/	
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            for(k=0; k<dimZ; k++) {
+				index = (dimX*dimY)*k + j*dimX+i;
+                /* boundary conditions */
+                i1 = i + 1; if (i == dimX-1) i1 = i;
+                j1 = j + 1; if (j == dimY-1) j1 = j;
+                k1 = k + 1; if (k == dimZ-1) k1 = k;
+                
+                /* Forward differences */                
+                NOMx_2 = powf((float)(U[(dimX*dimY)*k + j1*dimX+i] - U[index]),2); /* x+ */
+                NOMy_2 = powf((float)(U[(dimX*dimY)*k + j*dimX+i1] - U[index]),2); /* y+ */
+                NOMz_2 = powf((float)(U[(dimX*dimY)*k1 + j*dimX+i] - U[index]),2); /* z+ */
+                
+                E_Grad += sqrtf((float)(NOMx_2) + (float)(NOMy_2) + (float)(NOMz_2)); /* gradient term energy */
+                E_Data += 0.5f * lambda*(powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
+			}
+		}
+	}
+		E_val[0] = E_Grad + E_Data;
+		return *E_val;
+}
diff --git a/Core/regularisers_CPU/utils.h b/Core/regularisers_CPU/utils.h
index bd76bf0..fe08735 100644
--- a/Core/regularisers_CPU/utils.h
+++ b/Core/regularisers_CPU/utils.h
@@ -29,6 +29,7 @@ extern "C" {
 #endif
 CCPI_EXPORT float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
 CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY);
+CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY, int dimZ);
 #ifdef __cplusplus
 }
 #endif
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);
+    }
+}