summaryrefslogtreecommitdiffstats
path: root/main_func
diff options
context:
space:
mode:
Diffstat (limited to 'main_func')
-rw-r--r--main_func/FISTA_REC.m40
-rw-r--r--main_func/compile_mex.m11
-rw-r--r--main_func/regularizers_CPU/FGP_TV_core.c26
-rw-r--r--main_func/regularizers_CPU/FGP_TV_core.h27
-rw-r--r--main_func/regularizers_CPU/LLT_model.c26
-rw-r--r--main_func/regularizers_CPU/LLT_model_core.c25
-rw-r--r--main_func/regularizers_CPU/LLT_model_core.h25
-rw-r--r--main_func/regularizers_CPU/PatchBased_Regul.c43
-rw-r--r--main_func/regularizers_CPU/PatchBased_Regul_core.c57
-rw-r--r--main_func/regularizers_CPU/PatchBased_Regul_core.h33
-rw-r--r--main_func/regularizers_CPU/SplitBregman_TV_core.c1
-rw-r--r--main_func/regularizers_CPU/SplitBregman_TV_core.h24
-rw-r--r--main_func/regularizers_CPU/TGV_PD.c68
-rw-r--r--main_func/regularizers_CPU/TGV_PD_core.c1
-rw-r--r--main_func/regularizers_CPU/TGV_PD_core.h29
15 files changed, 180 insertions, 256 deletions
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 2823691..1e93719 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -99,7 +99,7 @@ else
clear proj_geomT vol_geomT
else
% divergen beam geometry
- fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry...');
+ fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry... will take some time!');
niter = 8; % number of iteration for PM
x1 = rand(N,N,SlicesZ);
sqweight = sqrt(weights);
@@ -216,6 +216,39 @@ if (isfield(params,'initialize'))
else
X = zeros(N,N,SlicesZ, 'single'); % storage for the solution
end
+if (isfield(params,'OS'))
+ % Ordered Subsets reorganisation of data and angles
+ OS = 1;
+ subsets = 8;
+ angles = proj_geom.ProjectionAngles;
+ binEdges = linspace(min(angles),max(angles),subsets+1);
+
+ % assign values to bins
+ [binsDiscr,~] = histc(angles, [binEdges(1:end-1) Inf]);
+
+ % rearrange angles into subsets
+ AnglesReorg = zeros(length(angles),1);
+ SinoReorg = zeros(Detectors, anglesNumb, SlicesZ, 'single');
+
+ counterM = 0;
+ for ii = 1:max(binsDiscr(:))
+ counter = 0;
+ for jj = 1:subsets
+ curr_index = ii+jj-1 + counter;
+ if (binsDiscr(jj) >= ii)
+ counterM = counterM + 1;
+ AnglesReorg(counterM) = angles(curr_index);
+ SinoReorg(:,counterM,:) = squeeze(sino(:,curr_index,:));
+ end
+ counter = (counter + binsDiscr(jj)) - 1;
+ end
+ end
+ sino = SinoReorg;
+ clear SinoReorg;
+else
+ OS = 0; % normal FISTA
+end
+
%----------------Reconstruction part------------------------
Resid_error = zeros(iterFISTA,1); % errors vector (if the ground truth is given)
@@ -241,7 +274,7 @@ for i = 1:iterFISTA
% ring removal part (Group-Huber fidelity)
for kkk = 1:anglesNumb
residual(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
- end
+ end
vec = sum(residual,2);
if (SlicesZ > 1)
vec = squeeze(vec(:,1,:));
@@ -249,7 +282,7 @@ for i = 1:iterFISTA
r = r_x - (1./L_const).*vec;
else
% no ring removal
- residual = weights.*(sino_updt - sino);
+ residual = weights.*(sino_updt - sino);
end
objective(i) = (0.5*norm(residual(:))^2)/(Detectors*anglesNumb*SlicesZ); % for the objective function output
@@ -305,5 +338,6 @@ end
output.Resid_error = Resid_error;
output.objective = objective;
output.L_const = L_const;
+output.sino = sino;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
diff --git a/main_func/compile_mex.m b/main_func/compile_mex.m
index 7bfa8eb..66c05da 100644
--- a/main_func/compile_mex.m
+++ b/main_func/compile_mex.m
@@ -1,10 +1,11 @@
-% compile mex's once
+% compile mex's in Matlab once
cd regularizers_CPU/
-mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex FGP_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex LLT_model.c LLT_model_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex FGP_TV.c FGP_TV_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex SplitBregman_TV.c SplitBregman_TV_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex TGV_PD.c TGV_PD_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex PatchBased_Regul.c PatchBased_Regul_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
cd ../../
cd demos
diff --git a/main_func/regularizers_CPU/FGP_TV_core.c b/main_func/regularizers_CPU/FGP_TV_core.c
index c383a71..a55991c 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.c
+++ b/main_func/regularizers_CPU/FGP_TV_core.c
@@ -19,6 +19,32 @@ limitations under the License.
#include "FGP_TV_core.h"
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Number of iterations [OPTIONAL parameter]
+ * 4. eplsilon: tolerance constant [OPTIONAL parameter]
+ * 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ * [2] last function value
+ *
+ * Example of image denoising:
+ * figure;
+ * Im = double(imread('lena_gray_256.tif'))/255; % loading image
+ * u0 = Im + .05*randn(size(Im)); % adding noise
+ * u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+ *
+ * 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"
+ *
+ * D. Kazantsev, 2016-17
+ *
+ */
+
/* 2D-case related Functions */
/*****************************************************************/
float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY)
diff --git a/main_func/regularizers_CPU/FGP_TV_core.h b/main_func/regularizers_CPU/FGP_TV_core.h
index 8d6af3e..e5328fb 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.h
+++ b/main_func/regularizers_CPU/FGP_TV_core.h
@@ -24,33 +24,6 @@ limitations under the License.
#include <stdio.h>
#include "omp.h"
-/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
-*
-* Input Parameters:
-* 1. Noisy image/volume [REQUIRED]
-* 2. lambda - regularization parameter [REQUIRED]
-* 3. Number of iterations [OPTIONAL parameter]
-* 4. eplsilon: tolerance constant [OPTIONAL parameter]
-* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
-*
-* Output:
-* [1] Filtered/regularized image
-* [2] last function value
-*
-* Example of image denoising:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255; % loading image
-* u0 = Im + .05*randn(size(Im)); % adding noise
-* u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-*
-* to compile with OMP support: mex FGP_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* 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"
-*
-* D. Kazantsev, 2016-17
-*
-*/
-
float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
diff --git a/main_func/regularizers_CPU/LLT_model.c b/main_func/regularizers_CPU/LLT_model.c
index 6b50d33..47146a5 100644
--- a/main_func/regularizers_CPU/LLT_model.c
+++ b/main_func/regularizers_CPU/LLT_model.c
@@ -20,6 +20,32 @@ limitations under the License.
#include "mex.h"
#include "LLT_model_core.h"
+/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
+*
+* Input Parameters:
+* 1. U0 - origanal noise image/volume
+* 2. lambda - regularization parameter
+* 3. tau - time-step for explicit scheme
+* 4. iter - iterations number
+* 5. epsil - tolerance constant (to terminate earlier)
+* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255; % loading image
+* u0 = Im + .03*randn(size(Im)); % adding noise
+* [Den] = LLT_model(single(u0), 10, 0.1, 1);
+*
+*
+* to compile with OMP support: mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
+*
+* 28.11.16/Harwell
+*/
+
void mexFunction(
int nlhs, mxArray *plhs[],
diff --git a/main_func/regularizers_CPU/LLT_model_core.c b/main_func/regularizers_CPU/LLT_model_core.c
index 3098269..7a1cdbe 100644
--- a/main_func/regularizers_CPU/LLT_model_core.c
+++ b/main_func/regularizers_CPU/LLT_model_core.c
@@ -19,6 +19,31 @@ limitations under the License.
#include "LLT_model_core.h"
+/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
+*
+* Input Parameters:
+* 1. U0 - origanal noise image/volume
+* 2. lambda - regularization parameter
+* 3. tau - time-step for explicit scheme
+* 4. iter - iterations number
+* 5. epsil - tolerance constant (to terminate earlier)
+* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255; % loading image
+* u0 = Im + .03*randn(size(Im)); % adding noise
+* [Den] = LLT_model(single(u0), 10, 0.1, 1);
+*
+* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
+*
+* 28.11.16/Harwell
+*/
+
+
float der2D(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ)
{
int i, j, i_p, i_m, j_m, j_p;
diff --git a/main_func/regularizers_CPU/LLT_model_core.h b/main_func/regularizers_CPU/LLT_model_core.h
index ef8803d..10f52fe 100644
--- a/main_func/regularizers_CPU/LLT_model_core.h
+++ b/main_func/regularizers_CPU/LLT_model_core.h
@@ -26,31 +26,6 @@ limitations under the License.
#define EPS 0.01
-/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
-*
-* Input Parameters:
-* 1. U0 - origanal noise image/volume
-* 2. lambda - regularization parameter
-* 3. tau - time-step for explicit scheme
-* 4. iter - iterations number
-* 5. epsil - tolerance constant (to terminate earlier)
-* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255; % loading image
-* u0 = Im + .03*randn(size(Im)); % adding noise
-* [Den] = LLT_model(single(u0), 10, 0.1, 1);
-*
-*
-* to compile with OMP support: mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
-*
-* 28.11.16/Harwell
-*/
/* 2D functions */
float der2D(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ);
float div_upd2D(float *U0, float *U, float *D1, float *D2, int dimX, int dimY, int dimZ, float lambda, float tau);
diff --git a/main_func/regularizers_CPU/PatchBased_Regul.c b/main_func/regularizers_CPU/PatchBased_Regul.c
index 24ee210..59eb3b3 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul.c
+++ b/main_func/regularizers_CPU/PatchBased_Regul.c
@@ -27,27 +27,23 @@ limitations under the License.
* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
* 2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
*
- * Input Parameters (mandatory):
- * 1. Image (2D or 3D)
- * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
- * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
- * 4. h - parameter for the PB penalty function
- * 5. lambda - regularization parameter
+ * Input Parameters:
+ * 1. Image (2D or 3D) [required]
+ * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window) [optional]
+ * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window) [optional]
+ * 4. h - parameter for the PB penalty function [optional]
+ * 5. lambda - regularization parameter [optional]
* Output:
* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
*
- * Quick 2D denoising example in Matlab:
+ * 2D denoising example in Matlab:
Im = double(imread('lena_gray_256.tif'))/255; % loading image
u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
- ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
- *
- * Please see more tests in a file:
- TestTemporalSmoothing.m
-
+ ImDen = PatchBased_Regul(single(u0), 3, 1, 0.08, 0.05);
*
* Matlab + C/mex compilers needed
- * to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
+ * to compile with OMP support: mex PatchBased_Regul.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
*
* D. Kazantsev *
* 02/07/2014
@@ -70,17 +66,23 @@ void mexFunction(
M = dims[1];
Z = dims[2];
- if ((numdims < 2) || (numdims > 3)) {mexErrMsgTxt("The input should be 2D image or 3D volume");}
+ if ((numdims < 2) || (numdims > 3)) {mexErrMsgTxt("The input is 2D image or 3D volume");}
if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input in single precision is required"); }
if(nrhs != 5) mexErrMsgTxt("Five inputs reqired: Image(2D,3D), SearchW, SimilW, Threshold, Regularization parameter");
/*Handling inputs*/
- A = (float *) mxGetData(prhs[0]); /* the image to regularize/filter */
- SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */
- SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */
- h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */
- lambda = (float) mxGetScalar(prhs[4]); /* regularization parameter */
+ A = (float *) mxGetData(prhs[0]); /* the image/volume to regularize/filter */
+ SearchW_real = 3; /*default value*/
+ SimilW = 1; /*default value*/
+ h = 0.1;
+ lambda = 0.1;
+
+ if ((nrhs == 2) || (nrhs == 3) || (nrhs == 4) || (nrhs == 5)) SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */
+ if ((nrhs == 4) || (nrhs == 5)) h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */
+ if ((nrhs == 5)) lambda = (float) mxGetScalar(prhs[4]); /* regularization parameter */
+
if (h <= 0) mexErrMsgTxt("Parmeter for the PB penalty function should be > 0");
if (lambda <= 0) mexErrMsgTxt(" Regularization parmeter should be > 0");
@@ -89,7 +91,6 @@ void mexFunction(
/* SearchW_full = 2*SearchW + 1; */ /* the full searching window size */
/* SimilW_full = 2*SimilW + 1; */ /* the full similarity window size */
-
padXY = SearchW + 2*SimilW; /* padding sizes */
newsizeX = N + 2*(padXY); /* the X size of the padded array */
@@ -135,4 +136,4 @@ void mexFunction(
switchpad_crop = 1; /*cropping*/
pad_crop(Bp, B, M, N, Z, newsizeY, newsizeX, newsizeZ, padXY, switchpad_crop);
} /*end else ndims*/
-} \ No newline at end of file
+}
diff --git a/main_func/regularizers_CPU/PatchBased_Regul_core.c b/main_func/regularizers_CPU/PatchBased_Regul_core.c
index 6f0a48d..acfb464 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul_core.c
+++ b/main_func/regularizers_CPU/PatchBased_Regul_core.c
@@ -19,40 +19,33 @@ limitations under the License.
#include "PatchBased_Regul_core.h"
-/* C-OMP implementation of patch-based (PB) regularization (2D and 3D cases).
-* This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
-*
-* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
-* 2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
-*
-* Input Parameters (mandatory):
-* 1. Image (2D or 3D)
-* 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
-* 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
-* 4. h - parameter for the PB penalty function
-* 5. lambda - regularization parameter
+/* C-OMP implementation of patch-based (PB) regularization (2D and 3D cases).
+ * This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
+ *
+ * References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
+ * 2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
+ *
+ * Input Parameters:
+ * 1. Image (2D or 3D) [required]
+ * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window) [optional]
+ * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window) [optional]
+ * 4. h - parameter for the PB penalty function [optional]
+ * 5. lambda - regularization parameter [optional]
-* Output:
-* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
-*
-* Quick 2D denoising example in Matlab:
-Im = double(imread('lena_gray_256.tif'))/255; % loading image
-u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-*
-* Please see more tests in a file:
-TestTemporalSmoothing.m
-
-*
-* Matlab + C/mex compilers needed
-* to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
-*
-* D. Kazantsev *
-* 02/07/2014
-* Harwell, UK
-*/
+ * Output:
+ * 1. regularized (denoised) Image (N x N)/volume (N x N x N)
+ *
+ * 2D denoising example in Matlab:
+ Im = double(imread('lena_gray_256.tif'))/255; % loading image
+ u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+ ImDen = PatchBased_Regul(single(u0), 3, 1, 0.08, 0.05);
+
+ * D. Kazantsev *
+ * 02/07/2014
+ * Harwell, UK
+ */
-/*2D version*/
+/*2D version function */
float PB_FUNC2D(float *A, float *B, int dimX, int dimY, int padXY, int SearchW, int SimilW, float h, float lambda)
{
int i, j, i_n, j_n, i_m, j_m, i_p, j_p, i_l, j_l, i1, j1, i2, j2, i3, j3, i5,j5, count, SimilW_full;
diff --git a/main_func/regularizers_CPU/PatchBased_Regul_core.h b/main_func/regularizers_CPU/PatchBased_Regul_core.h
index b83cf10..5aa6415 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul_core.h
+++ b/main_func/regularizers_CPU/PatchBased_Regul_core.h
@@ -26,39 +26,6 @@ limitations under the License.
#include <stdio.h>
#include "omp.h"
-/* C-OMP implementation of patch-based (PB) regularization (2D and 3D cases).
-* This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
-*
-* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
-* 2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
-*
-* Input Parameters (mandatory):
-* 1. Image (2D or 3D)
-* 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
-* 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
-* 4. h - parameter for the PB penalty function
-* 5. lambda - regularization parameter
-
-* Output:
-* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
-*
-* Quick 2D denoising example in Matlab:
-Im = double(imread('lena_gray_256.tif'))/255; % loading image
-u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-*
-* Please see more tests in a file:
-TestTemporalSmoothing.m
-
-*
-* Matlab + C/mex compilers needed
-* to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
-*
-* D. Kazantsev *
-* 02/07/2014
-* Harwell, UK
-*/
-
float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop);
float PB_FUNC2D(float *A, float *B, int dimX, int dimY, int padXY, int SearchW, int SimilW, float h, float lambda);
float PB_FUNC3D(float *A, float *B, int dimX, int dimY, int dimZ, int padXY, int SearchW, int SimilW, float h, float lambda);
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.c b/main_func/regularizers_CPU/SplitBregman_TV_core.c
index 26ad5b1..283dd43 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.c
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.c
@@ -37,7 +37,6 @@ limitations under the License.
* u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
* u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
*
-* to compile with OMP support: mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
* References:
* The Split Bregman Method for L1 Regularized Problems, by Tom Goldstein and Stanley Osher.
* D. Kazantsev, 2016*
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.h b/main_func/regularizers_CPU/SplitBregman_TV_core.h
index ed9112c..a7aaabb 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.h
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.h
@@ -23,30 +23,6 @@ limitations under the License.
#include <stdio.h>
#include "omp.h"
-/* C-OMP implementation of Split Bregman - TV denoising-regularization model (2D/3D)
-*
-* Input Parameters:
-* 1. Noisy image/volume
-* 2. lambda - regularization parameter
-* 3. Number of iterations [OPTIONAL parameter]
-* 4. eplsilon - tolerance constant [OPTIONAL parameter]
-* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255; % loading image
-* u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-* u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-*
-* to compile with OMP support: mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References:
-* The Split Bregman Method for L1 Regularized Problems, by Tom Goldstein and Stanley Osher.
-* D. Kazantsev, 2016*
-*/
-
float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
float gauss_seidel2D(float *U, float *A, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda, float mu);
float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda);
diff --git a/main_func/regularizers_CPU/TGV_PD.c b/main_func/regularizers_CPU/TGV_PD.c
index 8bce18a..d4bb787 100644
--- a/main_func/regularizers_CPU/TGV_PD.c
+++ b/main_func/regularizers_CPU/TGV_PD.c
@@ -39,7 +39,7 @@ limitations under the License.
* u0 = Im + .03*randn(size(Im)); % adding noise
* tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
*
- * to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+ * to compile with OMP support: mex TGV_PD.c TGV_PD_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
* References:
* K. Bredies "Total Generalized Variation"
*
@@ -53,7 +53,7 @@ void mexFunction(
{
int number_of_dims, iter, dimX, dimY, dimZ, ll;
const int *dim_array;
- float *A, *U, *U_old, *P1, *P2, *P3, *Q1, *Q2, *Q3, *Q4, *Q5, *Q6, *Q7, *Q8, *Q9, *V1, *V1_old, *V2, *V2_old, *V3, *V3_old, lambda, L2, tau, sigma, alpha1, alpha0;
+ float *A, *U, *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, lambda, L2, tau, sigma, alpha1, alpha0;
number_of_dims = mxGetNumberOfDimensions(prhs[0]);
dim_array = mxGetDimensions(prhs[0]);
@@ -88,40 +88,10 @@ void mexFunction(
V1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
V1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
V2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- V2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); }
- else if (number_of_dims == 3) {
- mexErrMsgTxt("The input data should be 2D");
- /*3D case*/
-// dimZ = dim_array[2];
-// U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//
-// P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//
-// Q1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q4 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q5 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q6 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q7 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q8 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// Q9 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//
-// U_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//
-// V1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// V1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// V2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// V2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// V3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-// V3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- }
- else {mexErrMsgTxt("The input data should be 2D");}
-
-
- /*printf("%i \n", i);*/
+ V2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+
+
+ /*printf("%i \n", i);*/
L2 = 12.0; /*Lipshitz constant*/
tau = 1.0/pow(L2,0.5);
sigma = 1.0/pow(L2,0.5);
@@ -129,7 +99,6 @@ void mexFunction(
/*Copy A to U*/
copyIm(A, U, dimX, dimY, dimZ);
- if (number_of_dims == 2) {
/* Here primal-dual iterations begin for 2D */
for(ll = 0; ll < iter; ll++) {
@@ -164,23 +133,12 @@ void mexFunction(
/*get new V*/
newU(V1, V1_old, dimX, dimY, dimZ);
newU(V2, V2_old, dimX, dimY, dimZ);
- } /*end of iterations*/
+ } /*end of iterations*/
}
-
-// /*3D version*/
-// if (number_of_dims == 3) {
-// /* Here primal-dual iterations begin for 3D */
-// for(ll = 0; ll < iter; ll++) {
-//
-// /* Calculate Dual Variable P */
-// DualP_3D(U, V1, V2, V3, P1, P2, P3, dimX, dimY, dimZ, sigma);
-//
-// /*Projection onto convex set for P*/
-// ProjP_3D(P1, P2, P3, dimX, dimY, dimZ, alpha1);
-//
-// /* Calculate Dual Variable Q */
-// DualQ_3D(V1, V2, V2, Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, Q9, dimX, dimY, dimZ, sigma);
-//
-// } /*end of iterations*/
-// }
+ else if (number_of_dims == 3) {
+ mexErrMsgTxt("The input data should be a 2D array");
+ /*3D case*/
+ }
+ else {mexErrMsgTxt("The input data should be a 2D array");}
+
}
diff --git a/main_func/regularizers_CPU/TGV_PD_core.c b/main_func/regularizers_CPU/TGV_PD_core.c
index 4c0427c..1164b73 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.c
+++ b/main_func/regularizers_CPU/TGV_PD_core.c
@@ -38,7 +38,6 @@ limitations under the License.
* u0 = Im + .03*randn(size(Im)); % adding noise
* tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
*
- * to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
* References:
* K. Bredies "Total Generalized Variation"
*
diff --git a/main_func/regularizers_CPU/TGV_PD_core.h b/main_func/regularizers_CPU/TGV_PD_core.h
index e25ae68..04ba95c 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.h
+++ b/main_func/regularizers_CPU/TGV_PD_core.h
@@ -24,32 +24,6 @@ limitations under the License.
#include <stdio.h>
#include "omp.h"
-/* C-OMP implementation of Primal-Dual denoising method for
-* Total Generilized Variation (TGV)-L2 model (2D case only)
-*
-* Input Parameters:
-* 1. Noisy image/volume (2D)
-* 2. lambda - regularization parameter
-* 3. parameter to control first-order term (alpha1)
-* 4. parameter to control the second-order term (alpha0)
-* 5. Number of CP iterations
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255; % loading image
-* u0 = Im + .03*randn(size(Im)); % adding noise
-* tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
-*
-* to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References:
-* K. Bredies "Total Generalized Variation"
-*
-* 28.11.16/Harwell
-*/
-
/* 2D functions */
float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, int dimZ, float sigma);
float ProjP_2D(float *P1, float *P2, int dimX, int dimY, int dimZ, float alpha1);
@@ -57,8 +31,5 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX,
float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, float alpha0);
float DivProjP_2D(float *U, float *A, float *P1, float *P2, int dimX, int dimY, int dimZ, float lambda, float tau);
float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, float tau);
-/*3D functions*/
-float DualP_3D(float *U, float *V1, float *V2, float *V3, float *P1, float *P2, float *P3, int dimX, int dimY, int dimZ, float sigma);
-
float newU(float *U, float *U_old, int dimX, int dimY, int dimZ);
float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 10:02:09 +0000