summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Core/regularisers_CPU/TGV_core.c202
-rw-r--r--Core/regularisers_CPU/TGV_core.h2
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m15
-rw-r--r--Wrappers/Matlab/demos/demoMatlab_denoise.m16
-rw-r--r--Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c29
5 files changed, 170 insertions, 94 deletions
diff --git a/Core/regularisers_CPU/TGV_core.c b/Core/regularisers_CPU/TGV_core.c
index 964c8a1..01306f3 100644
--- a/Core/regularisers_CPU/TGV_core.c
+++ b/Core/regularisers_CPU/TGV_core.c
@@ -121,7 +121,7 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in
DualP_3D(U, V1, V2, V3, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), sigma);
/*Projection onto convex set for P*/
- ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), alpha1);
+ ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), alpha1);
/* Calculate Dual Variable Q */
DualQ_3D(V1, V2, V3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), sigma);
@@ -189,10 +189,9 @@ float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1)
#pragma omp parallel for shared(P1,P2) private(i,j,index,grad_magn)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- grad_magn = sqrt(pow(P1[index],2) + pow(P2[index],2));
- grad_magn = grad_magn/alpha1;
- if (grad_magn > 1.0) {
+ index = j*dimX+i;
+ grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2)))/alpha1;
+ if (grad_magn > 1.0f) {
P1[index] /= grad_magn;
P2[index] /= grad_magn;
}
@@ -207,21 +206,14 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX,
#pragma omp parallel for shared(Q1,Q2,Q3,V1,V2) private(i,j,index,q1,q2,q11,q22)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- /* symmetric boundary conditions (Neuman) */
- if (i == dimX-1)
- { q1 = (V1[j*dimX+(i-1)] - V1[index]);
- q11 = (V2[j*dimX+(i-1)] - V2[index]);
- }
- else {
- q1 = (V1[j*dimX+(i+1)] - V1[index]);
- q11 = (V2[j*dimX+(i+1)] - V2[index]);
+ index = j*dimX+i;
+ q1 = 0.0f; q11 = 0.0f; q2 = 0.0f; q22 = 0.0f;
+ /* boundary conditions (Neuman) */
+ if (i != dimX-1){
+ q1 = V1[j*dimX+(i+1)] - V1[index];
+ q11 = V2[j*dimX+(i+1)] - V2[index];
}
- if (j == dimY-1) {
- q2 = (V2[(j-1)*dimX+i] - V2[index]);
- q22 = (V1[(j-1)*dimX+i] - V1[index]);
- }
- else {
+ if (j != dimY-1) {
q2 = V2[(j+1)*dimX+i] - V2[index];
q22 = V1[(j+1)*dimX+i] - V1[index];
}
@@ -238,10 +230,10 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alph
#pragma omp parallel for shared(Q1,Q2,Q3) private(i,j,index,grad_magn)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- grad_magn = sqrt(pow(Q1[index],2) + pow(Q2[index],2) + 2*pow(Q3[index],2));
+ index = j*dimX+i;
+ grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + 2*pow(Q3[index],2));
grad_magn = grad_magn/alpha0;
- if (grad_magn > 1.0) {
+ if (grad_magn > 1.0f) {
Q1[index] /= grad_magn;
Q2[index] /= grad_magn;
Q3[index] /= grad_magn;
@@ -279,37 +271,29 @@ float newU(float *U, float *U_old, long dimX, long dimY)
float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float tau)
{
long i, j, index;
- float q1, q11, q2, q22, div1, div2;
-#pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q11, q2, q22, div1, div2)
+ float q1, q3_x, q3_y, q2, div1, div2;
+#pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q3_x, q3_y, q2, div1, div2)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- /* symmetric boundary conditions (Neuman) */
- if (i == 0) {
- q1 = Q1[index];
- q11 = Q3[index];
- }
- else {
+ index = j*dimX+i;
+ q2 = 0.0f; q3_y = 0.0f; q1 = 0.0f; q3_x = 0.0;
+ /* boundary conditions (Neuman) */
+ if (i != 0) {
q1 = Q1[index] - Q1[j*dimX+(i-1)];
- q11 = Q3[index] - Q3[j*dimX+(i-1)];
+ q3_x = Q3[index] - Q3[j*dimX+(i-1)];
}
- if (j == 0) {
- q2 = Q2[index];
- q22 = Q3[index];
- }
- else {
+ if (j != 0) {
q2 = Q2[index] - Q2[(j-1)*dimX+i];
- q22 = Q3[index] - Q3[(j-1)*dimX+i];
+ q3_y = Q3[index] - Q3[(j-1)*dimX+i];
}
- div1 = q1 + q22;
- div2 = q2 + q11;
+ div1 = q1 + q3_y;
+ div2 = q3_x + q2;
V1[index] += tau*(P1[index] + div1);
V2[index] += tau*(P2[index] + div2);
}}
return 1;
}
-
/********************************************************************/
/***************************3D Functions*****************************/
/********************************************************************/
@@ -320,16 +304,14 @@ float DualP_3D(float *U, float *V1, float *V2, float *V3, float *P1, float *P2,
#pragma omp parallel for shared(U,V1,V2,V3,P1,P2,P3) private(i,j,k,index)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
- for(k=0; k<dimZ; k++) {
-
- index = (dimX*dimY)*k + j*dimX+i;
-
+ for(k=0; k<dimZ; k++) {
+ index = (dimX*dimY)*k + j*dimX+i;
/* symmetric boundary conditions (Neuman) */
if (i == dimX-1) P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i-1)] - U[index]) - V1[index]);
else P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index]) - V1[index]);
if (j == dimY-1) P2[index] += sigma*((U[(dimX*dimY)*k + (j-1)*dimX+i] - U[index]) - V2[index]);
else P2[index] += sigma*((U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index]) - V2[index]);
- if (j == dimZ-1) P3[index] += sigma*((U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index]) - V3[index]);
+ if (k == dimZ-1) P3[index] += sigma*((U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index]) - V3[index]);
else P3[index] += sigma*((U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index]) - V3[index]);
}}}
return 1;
@@ -344,8 +326,7 @@ float ProjP_3D(float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ,
for(j=0; j<dimY; j++) {
for(k=0; k<dimZ; k++) {
index = (dimX*dimY)*k + j*dimX+i;
- grad_magn = sqrt(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2));
- grad_magn = grad_magn/alpha1;
+ grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2)))/alpha1;
if (grad_magn > 1.0f) {
P1[index] /= grad_magn;
P2[index] /= grad_magn;
@@ -357,39 +338,29 @@ float ProjP_3D(float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ,
/*Calculating dual variable Q (using forward differences)*/
float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float sigma)
{
- long i,j,k, index;
- float q1, q2, q3, q4, q5, q6, q11, q22, q33, q44, q55, q66;
-#pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q4,q5,q6,q11,q22,q33,q44,q55,q66)
+ long i,j,k,index;
+ float q1, q2, q3, q11, q22, q33, q44, q55, q66;
+#pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q11,q22,q33,q44,q55,q66)
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
for(k=0; k<dimZ; k++) {
index = (dimX*dimY)*k + j*dimX+i;
+ q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
/* symmetric boundary conditions (Neuman) */
- if (i == dimX-1){
- q1 = V1[(dimX*dimY)*k + j*dimX+(i-1)] - V1[index];
- q11 = V2[(dimX*dimY)*k + j*dimX+(i-1)] - V2[index];
- q33 = V3[(dimX*dimY)*k + j*dimX+(i-1)] - V3[index];
- }
- else {
- q1 = (V1[(dimX*dimY)*k + j*dimX+(i+1)] - V1[index]);
- q11 = (V2[(dimX*dimY)*k + j*dimX+(i+1)] - V2[index]);
+ if (i != dimX-1){
+ q1 = V1[(dimX*dimY)*k + j*dimX+(i+1)] - V1[index];
+ q11 = V2[(dimX*dimY)*k + j*dimX+(i+1)] - V2[index];
q33 = V3[(dimX*dimY)*k + j*dimX+(i+1)] - V3[index];
}
- if (j == dimY-1) {
- q2 = (V2[(dimX*dimY)*k + (j-1)*dimX+i] - V2[index]);
- q22 = (V1[(dimX*dimY)*k + (j-1)*dimX+i] - V1[index]);
- }
- else {
+ if (j != dimY-1) {
q2 = V2[(dimX*dimY)*k + (j+1)*dimX+i] - V2[index];
q22 = V1[(dimX*dimY)*k + (j+1)*dimX+i] - V1[index];
+ q55 = V3[(dimX*dimY)*k + (j+1)*dimX+i] - V3[index];
}
- if (k == dimZ-1) {
- q3 = V3[(dimX*dimY)*(k-1) + j*dimX+i] - V3[index];
- q44 = V1[(dimX*dimY)*(k-1) + j*dimX+i] - V1[index];
- }
- else {
+ if (k != dimZ-1) {
q3 = V3[(dimX*dimY)*(k+1) + j*dimX+i] - V3[index];
q44 = V1[(dimX*dimY)*(k+1) + j*dimX+i] - V1[index];
+ q66 = V2[(dimX*dimY)*(k+1) + j*dimX+i] - V2[index];
}
Q1[index] += sigma*(q1);
@@ -397,8 +368,99 @@ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3,
Q3[index] += sigma*(q3);
Q4[index] += sigma*(0.5f*(q11 + q22));
Q5[index] += sigma*(0.5f*(q33 + q44));
+ Q6[index] += sigma*(0.5f*(q55 + q66));
}}}
return 1;
}
-
+float ProjQ_3D(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float alpha0)
+{
+ float grad_magn;
+ long i,j,k,index;
+#pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,grad_magn)
+ for(i=0; i<dimX; i++) {
+ for(j=0; j<dimY; j++) {
+ for(k=0; k<dimZ; k++) {
+ index = (dimX*dimY)*k + j*dimX+i;
+ grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + pow(Q3[index],2) + 2.0f*pow(Q4[index],2) + 2.0f*pow(Q5[index],2) + 2.0f*pow(Q6[index],2));
+ grad_magn = grad_magn/alpha0;
+ if (grad_magn > 1.0f) {
+ Q1[index] /= grad_magn;
+ Q2[index] /= grad_magn;
+ Q3[index] /= grad_magn;
+ Q4[index] /= grad_magn;
+ Q5[index] /= grad_magn;
+ Q6[index] /= grad_magn;
+ }
+ }}}
+ return 1;
+}
+/* Divergence and projection for P*/
+float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lambda, float tau)
+{
+ long i,j,k,index;
+ float P_v1, P_v2, P_v3, div;
+#pragma omp parallel for shared(U,U0,P1,P2,P3) private(i,j,k,index,P_v1,P_v2,P_v3,div)
+ for(i=0; i<dimX; i++) {
+ for(j=0; j<dimY; j++) {
+ for(k=0; k<dimZ; k++) {
+ index = (dimX*dimY)*k + j*dimX+i;
+ if (i == 0) P_v1 = P1[index];
+ else P_v1 = P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)];
+ if (j == 0) P_v2 = P2[index];
+ else P_v2 = P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i];
+ if (k == 0) P_v3 = P3[index];
+ else P_v3 = P3[index] - P3[(dimX*dimY)*(k-1) + (j)*dimX+i];
+
+ div = P_v1 + P_v2 + P_v3;
+ U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
+ }}}
+ return *U;
+}
+/*get update for V*/
+float UpdV_3D(float *V1, float *V2, float *V3, float *P1, float *P2, float *P3, float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float tau)
+{
+ long i,j,k,index;
+ float q1, q4x, q6x, q2, q4y, q5y, q6z, q5z, q3, div1, div2, div3;
+#pragma omp parallel for shared(V1,V2,V3,P1,P2,P3,Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,q1,q4x,q6x,q2,q4y,q5y,q6z,q5z,q3,div1,div2,div3)
+ for(i=0; i<dimX; i++) {
+ for(j=0; j<dimY; j++) {
+ for(k=0; k<dimZ; k++) {
+ index = (dimX*dimY)*k + j*dimX+i;
+ q1 = 0.0f; q4x= 0.0f; q6x= 0.0f; q2= 0.0f; q4y= 0.0f; q5y= 0.0f; q6z= 0.0f; q5z= 0.0f; q3= 0.0f;
+ /* symmetric boundary conditions (Neuman) */
+ if (i != 0) {
+ q1 = Q1[index] - Q1[(dimX*dimY)*k + j*dimX+(i-1)];
+ q4x = Q4[index] - Q4[(dimX*dimY)*k + j*dimX+(i-1)];
+ q6x = Q6[index] - Q6[(dimX*dimY)*k + j*dimX+(i-1)];
+ }
+ if (j != 0) {
+ q2 = Q2[index] - Q2[(dimX*dimY)*k + (j-1)*dimX+i];
+ q4y = Q4[index] - Q4[(dimX*dimY)*k + (j-1)*dimX+i];
+ q5y = Q5[index] - Q5[(dimX*dimY)*k + (j-1)*dimX+i];
+ }
+ if (k != 0) {
+ q6z = Q6[index] - Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
+ q5z = Q5[index] - Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
+ q3 = Q3[index] - Q3[(dimX*dimY)*(k-1) + (j)*dimX+i];
+ }
+ div1 = q1 + q4y + q6z;
+ div2 = q4x + q2 + q5z;
+ div3 = q6x + q5y + q3;
+
+ V1[index] += tau*(P1[index] + div1);
+ V2[index] += tau*(P2[index] + div2);
+ V3[index] += tau*(P3[index] + div3);
+ }}}
+ return 1;
+}
+/*get updated solution U*/
+float newU3D(float *U, float *U_old, long dimX, long dimY, long dimZ)
+{
+ long i;
+#pragma omp parallel for shared(U,U_old) private(i)
+ for(i=0; i<dimX*dimY*dimZ; i++) {
+ U[i] = 2.0f*U[i] - U_old[i];
+ }
+ return *U;
+}
diff --git a/Core/regularisers_CPU/TGV_core.h b/Core/regularisers_CPU/TGV_core.h
index 8bd8cf6..8ff084e 100644
--- a/Core/regularisers_CPU/TGV_core.h
+++ b/Core/regularisers_CPU/TGV_core.h
@@ -65,7 +65,7 @@ CCPI_EXPORT float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2
CCPI_EXPORT float ProjQ_3D(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float alpha0);
CCPI_EXPORT float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lambda, float tau);
CCPI_EXPORT float UpdV_3D(float *V1, float *V2, float *V3, float *P1, float *P2, float *P3, float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float tau);
-CCPI_EXPORT float newU3D(float *U, float *U_old, long dimX, long dimY long dimZ);
+CCPI_EXPORT float newU3D(float *U, float *U_old, long dimX, long dimY, long dimZ);
#ifdef __cplusplus
}
#endif
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
index 5cc47b3..23cda32 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -2,11 +2,13 @@
clear; close all
Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+Path3 = sprintf(['..' filesep 'supp'], 1i);
addpath(Path1);
addpath(Path2);
+addpath(Path3);
N = 512;
-slices = 15;
+slices = 7;
vol3D = zeros(N,N,slices, 'single');
Ideal3D = zeros(N,N,slices, 'single');
Im = double(imread('lena_gray_512.tif'))/255; % loading image
@@ -131,7 +133,16 @@ figure; imshow(u_diff4(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (CP
% fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4);
% figure; imshow(u_diff4_g(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (GPU)');
%%
-
+fprintf('Denoise using the TGV model (CPU) \n');
+lambda_TGV = 0.05; % regularisation parameter
+alpha1 = 1.0; % parameter to control the first-order term
+alpha0 = 2.0; % parameter to control the second-order term
+iter_TGV = 40; % number of Primal-Dual iterations for TGV
+tic; u_tgv = TGV(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV, 128); toc;
+rmseTGV = RMSE(Ideal3D(:),u_tgv(:));
+fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV);
+figure; imshow(u_tgv(:,:,3), [0 1]); title('TGV denoised volume (CPU)');
+%%
%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< %
fprintf('Denoise a volume using the FGP-dTV model (CPU) \n');
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 3506cca..14d3096 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -60,20 +60,20 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)');
% figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)');
%%
fprintf('Denoise using the TGV model (CPU) \n');
-lambda_TGV = 0.04; % regularisation parameter
-alpha1 = 1; % parameter to control the first-order term
-alpha0 = 0.7; % parameter to control the second-order term
-iter_TGV = 500; % number of Primal-Dual iterations for TGV
+lambda_TGV = 0.045; % regularisation parameter
+alpha1 = 1.0; % parameter to control the first-order term
+alpha0 = 2.0; % parameter to control the second-order term
+iter_TGV = 2000; % number of Primal-Dual iterations for TGV
tic; u_tgv = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc;
rmseTGV = (RMSE(u_tgv(:),Im(:)));
fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV);
figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)');
%%
% fprintf('Denoise using the TGV model (GPU) \n');
-% lambda_TGV = 0.04; % regularisation parameter
-% alpha1 = 1; % parameter to control the first-order term
-% alpha0 = 0.7; % parameter to control the second-order term
-% iter_TGV = 500; % number of Primal-Dual iterations for TGV
+% lambda_TGV = 0.045; % regularisation parameter
+% alpha1 = 1.0; % parameter to control the first-order term
+% alpha0 = 2.0; % parameter to control the second-order term
+% iter_TGV = 2000; % number of Primal-Dual iterations for TGV
% tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc;
% rmseTGV_gpu = (RMSE(u_tgv_gpu(:),Im(:)));
% fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV_gpu);
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c
index 5459bf5..aa4eed4 100644
--- a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c
@@ -21,14 +21,14 @@ limitations under the License.
#include "TGV_core.h"
/* C-OMP implementation of Primal-Dual denoising method for
- * Total Generilized Variation (TGV)-L2 model [1] (2D case only)
+ * Total Generilized Variation (TGV)-L2 model [1] (2D/3D)
*
* Input Parameters:
- * 1. Noisy image (2D) (required)
- * 2. lambda - regularisation parameter (required)
- * 3. parameter to control the first-order term (alpha1) (default - 1)
- * 4. parameter to control the second-order term (alpha0) (default - 0.5)
- * 5. Number of Chambolle-Pock (Primal-Dual) iterations (default is 300)
+ * 1. Noisy image/volume (2D/3D)
+ * 2. lambda - regularisation parameter
+ * 3. parameter to control the first-order term (alpha1)
+ * 4. parameter to control the second-order term (alpha0)
+ * 5. Number of Chambolle-Pock (Primal-Dual) iterations
* 6. Lipshitz constant (default is 12)
*
* Output:
@@ -44,7 +44,7 @@ void mexFunction(
{
int number_of_dims, iter;
- mwSize dimX, dimY;
+ mwSize dimX, dimY, dimZ;
const mwSize *dim_array;
float *Input, *Output=NULL, lambda, alpha0, alpha1, L2;
@@ -55,7 +55,7 @@ void mexFunction(
/*Handling Matlab input data*/
if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant");
- Input = (float *) mxGetData(prhs[0]); /*noisy image (2D) */
+ Input = (float *) mxGetData(prhs[0]); /*noisy image/volume */
lambda = (float) mxGetScalar(prhs[1]); /* regularisation parameter */
alpha1 = 1.0f; /* parameter to control the first-order term */
alpha0 = 0.5f; /* parameter to control the second-order term */
@@ -69,12 +69,15 @@ void mexFunction(
if (nrhs == 6) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */
/*Handling Matlab output data*/
- dimX = dim_array[0]; dimY = dim_array[1];
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
if (number_of_dims == 2) {
- Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- /* running the function */
- TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY);
+ dimZ = 1; /*2D case*/
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
}
- if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");}
+ if (number_of_dims == 3) {
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ /* running the function */
+ TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-06 17:42:05 +0000