some corrections to energy estimation

5 files changed, 26 insertions, 22 deletions

diff --git a/Core/regularisers_CPU/utils.c b/Core/regularisers_CPU/utils.c
index f21d383..ca5c59a 100644
--- a/Core/regularisers_CPU/utils.c
+++ b/Core/regularisers_CPU/utils.c
@@ -36,10 +36,11 @@ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
 }
 */
 
-/* function that calculates TV energy (ROF model) 
- * min||\nabla u|| + 0.5*lambda*||u -u0||^2
+/* function that calculates TV energy
+ * type - 1:  2*lambda*min||\nabla u|| + ||u -u0||^2
+ * type - 2:  2*lambda*min||\nabla u|| 
  * */
-float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY)
+float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY)
 {
 	int i, j, i1, j1, index;
 	float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
@@ -55,15 +56,16 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int
                 /* Forward differences */                
                 NOMx_2 = powf((float)(U[j1*dimX + i] - U[index]),2); /* x+ */
                 NOMy_2 = powf((float)(U[j*dimX + i1] - U[index]),2); /* y+ */
-                E_Grad += sqrtf((float)(NOMx_2) + (float)(NOMy_2)); /* gradient term energy */
-                E_Data += 0.5f * lambda*(powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
+                E_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2)); /* gradient term energy */
+                E_Data += powf((float)(U[index]-U0[index]),2); /* fidelity term energy */
 			}
 		}
-		E_val[0] = E_Grad + E_Data;
+		if (type == 1) E_val[0] = E_Grad + E_Data;
+		if (type == 2) E_val[0] = E_Grad;
 		return *E_val;
 }
 
-float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY, int dimZ)
+float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, 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;
@@ -83,11 +85,12 @@ float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int dimX, int
                 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_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2) + (float)(NOMz_2)); /* gradient term energy */
+                E_Data += (powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
 			}
 		}
 	}
-		E_val[0] = E_Grad + E_Data;
+		if (type == 1) E_val[0] = E_Grad + E_Data;
+		if (type == 2) E_val[0] = E_Grad;
 		return *E_val;
 }
diff --git a/Core/regularisers_CPU/utils.h b/Core/regularisers_CPU/utils.h
index fe08735..866bc01 100644
--- a/Core/regularisers_CPU/utils.h
+++ b/Core/regularisers_CPU/utils.h
@@ -28,8 +28,8 @@ limitations under the License.
 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);
+CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY);
+CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, 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 84889d7..5a54d18 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -21,7 +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
+energyfunc_val_rof = TV_energy(single(u_rof),single(vol3D),lambda_reg, 1);  % 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');
@@ -34,7 +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
+energyfunc_val_fgp = TV_energy(single(u_fgp),single(vol3D),lambda_reg, 1); % 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');
@@ -47,7 +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
+energyfunc_val_sb = TV_energy(single(u_sb),single(vol3D),lambda_reg, 1);  % 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 526d21c..151a604 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -14,7 +14,7 @@ fprintf('Denoise using the ROF-TV model (CPU) \n');
 tau_rof = 0.0025; % time-marching constant 
 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
+energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg, 1);  % 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');
@@ -27,7 +27,7 @@ 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
+energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value
 figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)');
 
 %%
@@ -41,7 +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
+energyfunc_val_sb = TV_energy(single(u_sb),single(u0),lambda_reg, 1);  % 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/regularisers_CPU/TV_energy.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c
index 421bd4c..f9eb2ce 100644
--- a/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TV_energy.c
@@ -36,7 +36,7 @@ void mexFunction(
         int nrhs, const mxArray *prhs[])
         
 {
-    int number_of_dims, dimX, dimY, dimZ;
+    int number_of_dims, dimX, dimY, dimZ, type;
     const int  *dim_array;
     float *Input, *Input0, lambda;
     
@@ -44,11 +44,12 @@ void mexFunction(
     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");
+    if ((nrhs != 4)) mexErrMsgTxt("4 inputs: Two images or volumes of the same size required, estimated and the original (noisy), regularisation parameter, type");
     
     Input  = (float *) mxGetData(prhs[0]); /* Denoised Image/volume */
     Input0  = (float *) mxGetData(prhs[1]); /* Original (noisy) Image/volume */
     lambda =  (float) mxGetScalar(prhs[2]); /* regularisation parameter */
+    type =  (int) mxGetScalar(prhs[3]); /* type of energy */
     
     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"); }
@@ -61,9 +62,9 @@ void mexFunction(
     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);
+		TV_energy2D(Input, Input0, funcvalA, lambda, type, dimX, dimY);
 		}
     if (number_of_dims == 3) {
-        TV_energy3D(Input, Input0, funcvalA, lambda, dimX, dimY, dimZ);
+        TV_energy3D(Input, Input0, funcvalA, lambda, type, dimX, dimY, dimZ);
     }
 }