From a9d773b384c6391dbb9913deeafa3e79e108b790 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Fri, 20 Apr 2018 12:38:38 +0100
Subject: some corrections to energy estimation
---
Core/regularisers_CPU/utils.c | 23 ++++++++++++----------
Core/regularisers_CPU/utils.h | 4 ++--
Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 6 +++---
Wrappers/Matlab/demos/demoMatlab_denoise.m | 6 +++---
.../mex_compile/regularisers_CPU/TV_energy.c | 9 +++++----
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);
}
}
--
cgit v1.2.3