summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.c88
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.h7
-rw-r--r--src/Matlab/mex_compile/compileCPU_mex_Linux.m13
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c28
-rw-r--r--src/Python/ccpi/filters/regularisers.py7
-rw-r--r--src/Python/setup-regularisers.py.in2
-rw-r--r--src/Python/src/cpu_regularisers.pyx38
7 files changed, 102 insertions, 81 deletions
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index 68d58b7..f613f0e 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -28,28 +28,33 @@ limitations under the License.
* 4. eplsilon: tolerance constant
* 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
* 6. nonneg: 'nonnegativity (0 is OFF by default)
- * 7. print information: 0 (off) or 1 (on)
*
* Output:
- * [1] Filtered/regularized image
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
*
* This function is based on the Matlab's code and paper by
* [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
*/
-float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ)
+float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
{
- int ll;
+ int ll;
long j, DimTotal;
- float re, re1;
- float tk = 1.0f;
- float tkp1=1.0f;
+ float re, re1;
+ re = 0.0f; re1 = 0.0f;
+ float tk = 1.0f;
+ float tkp1 =1.0f;
int count = 0;
+
+ /*adding info into info_vector */
+// infovector[0] = 50.1f; /*iterations number (if stopped earlier based on tolerance)*/
+// infovector[1] = 0.55f; /* reached tolerance */
if (dimZ <= 1) {
- /*2D case */
- float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
- DimTotal = (long)(dimX*dimY);
+ /*2D case */
+ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
+ DimTotal = (long)(dimX*dimY);
Output_prev = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
@@ -59,7 +64,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
R1 = calloc(DimTotal, sizeof(float));
R2 = calloc(DimTotal, sizeof(float));
- /* begin iterations */
+ /* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
/* computing the gradient of the objective function */
@@ -79,24 +84,32 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
/* check early stopping criteria */
- re = 0.0f; re1 = 0.0f;
- for(j=0; j<DimTotal; j++)
- {
- re += pow(Output[j] - Output_prev[j],2);
- re1 += pow(Output[j],2);
+ if (epsil != 0.0f) {
+ for(j=0; j<DimTotal; j++)
+ {
+ re += pow(Output[j] - Output_prev[j],2);
+ re1 += pow(Output[j],2);
+ }
+ re = sqrt(re)/sqrt(re1);
+ if (re < epsil) count++;
+ if (count > 4) break;
+
+ copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
}
- re = sqrt(re)/sqrt(re1);
- if (re < epsil) count++;
- if (count > 4) break;
-
/*storing old values*/
- copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
tk = tkp1;
}
- if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll);
- free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
+ /*adding info into info_vector */
+ //info_vector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
+// info_vector[1] = 0.55f; /* reached tolerance */
+
+ copyIm(Input, infovector, (long)(dimX), (long)(dimY), 1l);
+
+// printf("%f\n", infovector[128]);
+
+ free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
}
else {
/*3D case*/
@@ -134,26 +147,31 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
/* calculate norm - stopping rules*/
- re = 0.0f; re1 = 0.0f;
- for(j=0; j<DimTotal; j++)
- {
- re += pow(Output[j] - Output_prev[j],2);
- re1 += pow(Output[j],2);
- }
+ if (epsil != 0.0f) {
+ for(j=0; j<DimTotal; j++)
+ {
+ re += pow(Output[j] - Output_prev[j],2);
+ re1 += pow(Output[j],2);
+ }
re = sqrt(re)/sqrt(re1);
/* stop if the norm residual is less than the tolerance EPS */
if (re < epsil) count++;
- if (count > 4) break;
-
- /*storing old values*/
- copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+ if (count > 4) break;
+
+ copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+ }
+
+ /*storing old values*/
copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
tk = tkp1;
}
- if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll);
- free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
+ /*adding info into info_vector */
+ //infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
+ //infovector[1] = re; /* reached tolerance */
+
+ free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
}
return *Output;
}
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h
index 3418604..04e6e80 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -35,10 +35,11 @@ limitations under the License.
* 4. eplsilon: tolerance constant
* 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
* 6. nonneg: 'nonnegativity (0 is OFF by default)
- * 7. print information: 0 (off) or 1 (on)
*
* Output:
- * [1] Filtered/regularized image
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+
*
* This function is based on the Matlab's code and paper by
* [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
@@ -47,7 +48,7 @@ limitations under the License.
#ifdef __cplusplus
extern "C" {
#endif
-CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 72a828e..f3d9ce1 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -2,9 +2,9 @@
fsep = '/';
-pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
-pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
-pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
+pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
+pathcopyFrom1 = sprintf(['..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
copyfile(pathcopyFrom, 'regularisers_CPU');
copyfile(pathcopyFrom1, 'regularisers_CPU');
@@ -76,6 +76,7 @@ delete PatchSelect_core* Nonlocal_TV_core*
delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
-pathA2 = sprintf(['..' fsep '..' fsep], 1i);
-cd(pathA2);
-cd demos
+%pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+% cd(pathA2);
+cd('/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/demos')
+%cd demos
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
index 642362f..603e0f4 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
@@ -47,13 +47,13 @@ void mexFunction(
int number_of_dims, iter, methTV, printswitch, nonneg;
mwSize dimX, dimY, dimZ;
const mwSize *dim_array;
- float *Input, *Output=NULL, lambda, epsil;
+ float *Input, *Output=NULL, *Output2=NULL, lambda, epsil;
number_of_dims = mxGetNumberOfDimensions(prhs[0]);
dim_array = mxGetDimensions(prhs[0]);
/*Handling Matlab input data*/
- if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+ if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch");
Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
@@ -61,27 +61,24 @@ void mexFunction(
epsil = 0.0001; /* default tolerance constant */
methTV = 0; /* default isotropic TV penalty */
nonneg = 0; /* default nonnegativity switch, off - 0 */
- printswitch = 0; /*default print is switched, off - 0 */
+// printswitch = 0; /*default print is switched, off - 0 */
if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
- if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
- if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
- if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) {
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6)) {
char *penalty_type;
penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
mxFree(penalty_type);
}
- if ((nrhs == 6) || (nrhs == 7)) {
+ if (nrhs == 6) {
nonneg = (int) mxGetScalar(prhs[5]);
if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
}
- if (nrhs == 7) {
- printswitch = (int) mxGetScalar(prhs[6]);
- if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
- }
+
/*Handling Matlab output data*/
dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
@@ -89,9 +86,12 @@ void mexFunction(
if (number_of_dims == 2) {
dimZ = 1; /*2D case*/
Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ Output2 = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
}
- if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
/* running the function */
- TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
-} \ No newline at end of file
+ TV_FGP_CPU_main(Input, Output, Output2, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ);
+}
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 588ea32..fb2c999 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -29,15 +29,14 @@ def ROF_TV(inputData, regularisation_parameter, iterations,
.format(device))
def FGP_TV(inputData, regularisation_parameter,iterations,
- tolerance_param, methodTV, nonneg, printM, device='cpu'):
+ tolerance_param, methodTV, nonneg, device='cpu'):
if device == 'cpu':
return TV_FGP_CPU(inputData,
regularisation_parameter,
iterations,
tolerance_param,
methodTV,
- nonneg,
- printM)
+ nonneg)
elif device == 'gpu' and gpu_enabled:
return TV_FGP_GPU(inputData,
regularisation_parameter,
@@ -45,7 +44,7 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
tolerance_param,
methodTV,
nonneg,
- printM)
+ 1)
else:
if not gpu_enabled and device == 'gpu':
raise ValueError ('GPU is not available')
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 82d9f9f..39b820a 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -44,7 +44,7 @@ extra_include_dirs += [os.path.join(".." , "Core"),
os.path.join(".." , "Core", "regularisers_GPU" , "LLTROF" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "NDF" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "dTV_FGP" ) ,
- os.path.join(".." , "Core", "regularisers_GPU" , "DIFF4th" ) ,
+ os.path.join(".." , "Core", "regularisers_GPU" , "Diff4th" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "PatchSelect" ) ,
"."]
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 11a0617..b7d029d 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -19,7 +19,7 @@ import numpy as np
cimport numpy as np
cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
-cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
cdef extern float SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ);
cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
@@ -45,7 +45,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste
def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularisation_parameter,
- int iterationsNumb,
+ int iterationsNumb,
float marching_step_parameter):
cdef long dims[2]
dims[0] = inputData.shape[0]
@@ -80,45 +80,46 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
#********************** Total-variation FGP *********************#
#****************************************************************#
#******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
+def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg):
if inputData.ndim == 2:
- return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+ return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg)
elif inputData.ndim == 3:
- return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+ return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg)
def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularisation_parameter,
int iterationsNumb,
float tolerance_param,
int methodTV,
- int nonneg,
- int printM):
-
+ int nonneg):
+
cdef long dims[2]
dims[0] = inputData.shape[0]
dims[1] = inputData.shape[1]
cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
np.zeros([dims[0],dims[1]], dtype='float32')
-
+
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] infovec = \
+ np.ones([dims[0],dims[1]], dtype='float32')
+
#/* Run FGP-TV iterations for 2D data */
- TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter,
+ TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0,0], regularisation_parameter,
iterationsNumb,
tolerance_param,
methodTV,
nonneg,
- printM,
dims[1],dims[0],1)
- return outputData
+ return (outputData,infovec)
def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
float regularisation_parameter,
int iterationsNumb,
float tolerance_param,
int methodTV,
- int nonneg,
- int printM):
+ int nonneg):
+
cdef long dims[3]
dims[0] = inputData.shape[0]
dims[1] = inputData.shape[1]
@@ -126,16 +127,17 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-
+ cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+ np.zeros([2], dtype='float32')
+
#/* Run FGP-TV iterations for 3D data */
- TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
iterationsNumb,
tolerance_param,
methodTV,
nonneg,
- printM,
dims[2], dims[1], dims[0])
- return outputData
+ return (outputData,infovec)
#***************************************************************#
#********************** Total-variation SB *********************#
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 23:49:50 +0000