From bc286e76483c366ba221ce79349a277fb6db32ed Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Fri, 26 Jan 2018 12:11:07 +0000
Subject: ROF TV regularizer added #22
---
Core/regularizers_CPU/ROF_TV_core.c | 259 +++++++++++++++++++++
Core/regularizers_CPU/ROF_TV_core.h | 54 +++++
Wrappers/Matlab/mex_compile/compile_mex.m | 3 +-
.../Matlab/mex_compile/regularizers_CPU/ROF_TV.c | 106 +++++++++
4 files changed, 421 insertions(+), 1 deletion(-)
create mode 100644 Core/regularizers_CPU/ROF_TV_core.c
create mode 100644 Core/regularizers_CPU/ROF_TV_core.h
create mode 100644 Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
diff --git a/Core/regularizers_CPU/ROF_TV_core.c b/Core/regularizers_CPU/ROF_TV_core.c
new file mode 100644
index 0000000..0b24806
--- /dev/null
+++ b/Core/regularizers_CPU/ROF_TV_core.c
@@ -0,0 +1,259 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "ROF_TV_core.h"
+
+#define EPS 0.000001
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+int sign(float x) {
+ return (x > 0) - (x < 0);
+}
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+
+/* calculate differences 1 */
+float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ)
+{
+ float NOMx_1, NOMy_1, NOMy_0, NOMz_1, NOMz_0, denom1, denom2,denom3, T1;
+ int i,j,k,i1,i2,k1,j1,j2,k2;
+
+ if (dimZ == 0) {
+#pragma omp parallel for shared (A, D1, dimX, dimY) private(i, j, i1, j1, i2, j2,NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+
+ /* Forward-backward differences */
+ NOMx_1 = A[i1*dimY + j] - A[(i)*dimY + j]; /* x+ */
+ NOMy_1 = A[i*dimY + j1] - A[(i)*dimY + j]; /* y+ */
+ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */
+ NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; /* y- */
+
+ denom1 = NOMx_1*NOMx_1;
+ denom2 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+ denom2 = denom2*denom2;
+ T1 = sqrt(denom1 + denom2 + EPS);
+ D1[i*dimY+j] = NOMx_1/T1;
+ }}
+ }
+ else {
+#pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1,NOMy_1,NOMy_0,NOMz_1,NOMz_0,denom1,denom2,denom3,T1)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ for(k=0; k<dimZ; k++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+ k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
+ k2 = k - 1; if (k2 < 0) k2 = k+1;
+ /*B[(dimX*dimY)*k + i*dimY+j] = 0.25*(A[(dimX*dimY)*k + (i1)*dimY + j] + A[(dimX*dimY)*k + (i2)*dimY + j] + A[(dimX*dimY)*k + (i)*dimY + j1] + A[(dimX*dimY)*k + (i)*dimY + j2]) - A[(dimX*dimY)*k + i*dimY + j];*/
+
+ /* Forward-backward differences */
+ NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+ NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */
+ NOMy_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j2]; /* y- */
+
+ NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* z+ */
+ NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; /* z- */
+
+
+ denom1 = NOMx_1*NOMx_1;
+ denom2 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+ denom2 = denom2*denom2;
+ denom3 = 0.5*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0)));
+ denom3 = denom3*denom3;
+ T1 = sqrt(denom1 + denom2 + denom3 + EPS);
+ D1[(dimX*dimY)*k + i*dimY+j] = NOMx_1/T1;
+ }}}
+ }
+ return *D1;
+}
+/* calculate differences 2 */
+float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ)
+{
+ float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2;
+ int i,j,k,i1,i2,k1,j1,j2,k2;
+
+ if (dimZ == 0) {
+#pragma omp parallel for shared (A, D2, dimX, dimY) private(i, j, i1, j1, i2, j2, NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+
+ /* Forward-backward differences */
+ NOMx_1 = A[(i1)*dimY + j] - A[(i)*dimY + j]; /* x+ */
+ NOMy_1 = A[i*dimY + j1] - A[(i)*dimY + j]; /* y+ */
+ NOMx_0 = A[(i)*dimY + j] - A[(i2)*dimY + j]; /* x- */
+ /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */ /* y- */
+
+ denom1 = NOMy_1*NOMy_1;
+ denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+ denom2 = denom2*denom2;
+ T2 = sqrt(denom1 + denom2 + EPS);
+ D2[i*dimY+j] = NOMy_1/T2;
+ }}
+ }
+ else {
+#pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ for(k=0; k<dimZ; k++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+ k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
+ k2 = k - 1; if (k2 < 0) k2 = k+1;
+
+
+ /* Forward-backward differences */
+ NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+ NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+ NOMx_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i2)*dimY + j]; /* x- */
+ NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* z+ */
+ NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; /* z- */
+
+
+ denom1 = NOMy_1*NOMy_1;
+ denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+ denom2 = denom2*denom2;
+ denom3 = 0.5*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0)));
+ denom3 = denom3*denom3;
+ T2 = sqrt(denom1 + denom2 + denom3 + EPS);
+ D2[(dimX*dimY)*k + i*dimY+j] = NOMy_1/T2;
+ }}}
+ }
+ return *D2;
+}
+
+/* calculate differences 3 */
+float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ)
+{
+ float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3;
+ int i,j,k,i1,i2,k1,j1,j2,k2;
+
+#pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMy_0, NOMx_0, NOMz_1, denom1, denom2, denom3, T3)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ for(k=0; k<dimZ; k++) {
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+ k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
+ k2 = k - 1; if (k2 < 0) k2 = k+1;
+
+ /* Forward-backward differences */
+ NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+ NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+ NOMy_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j2]; /* y- */
+ NOMx_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i2)*dimY + j]; /* x- */
+ NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* z+ */
+ /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */
+
+ denom1 = NOMz_1*NOMz_1;
+ denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+ denom2 = denom2*denom2;
+ denom3 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+ denom3 = denom3*denom3;
+ T3 = sqrt(denom1 + denom2 + denom3 + EPS);
+ D3[(dimX*dimY)*k + i*dimY+j] = NOMz_1/T3;
+ }}}
+ return *D3;
+}
+
+/* calculate divergence */
+float TV_main(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ)
+{
+ float dv1, dv2, dv3;
+ int index,i,j,k,i1,i2,k1,j1,j2,k2;
+
+ if (dimZ == 0) {
+#pragma omp parallel for shared (D1, D2, B, dimX, dimY) private(index, i, j, i1, j1, i2, j2,dv1,dv2)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ index = (i)*dimY + j;
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+
+ /* divergence components */
+ dv1 = D1[index] - D1[(i2)*dimY + j];
+ dv2 = D2[index] - D2[(i)*dimY + j2];
+
+ B[index] = B[index] + tau*lambda*(dv1 + dv2) + tau*(A[index] - B[index]);
+
+ }}
+ }
+ else {
+#pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ for(k=0; k<dimZ; k++) {
+ index = (dimX*dimY)*k + i*dimY+j;
+ /* symmetric boundary conditions (Neuman) */
+ i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+ i2 = i - 1; if (i2 < 0) i2 = i+1;
+ j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+ j2 = j - 1; if (j2 < 0) j2 = j+1;
+ k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
+ k2 = k - 1; if (k2 < 0) k2 = k+1;
+
+ /*divergence components */
+ dv1 = D1[index] - D1[(dimX*dimY)*k + i2*dimY+j];
+ dv2 = D2[index] - D2[(dimX*dimY)*k + i*dimY+j2];
+ dv3 = D3[index] - D3[(dimX*dimY)*k2 + i*dimY+j];
+
+ B[index] = B[index] + tau*lambda*(dv1 + dv2 + dv3) + tau*(A[index] - B[index]);
+ }}}
+ }
+ return *B;
+}
\ No newline at end of file
diff --git a/Core/regularizers_CPU/ROF_TV_core.h b/Core/regularizers_CPU/ROF_TV_core.h
new file mode 100644
index 0000000..6e4f961
--- /dev/null
+++ b/Core/regularizers_CPU/ROF_TV_core.h
@@ -0,0 +1,54 @@
+/*
+This work is part of the Core Imaging Library developed by
+Visual Analytics and Imaging System Group of the Science Technology
+Facilities Council, STFC
+
+Copyright 2017 Daniil Kazantsev
+Copyright 2017 Srikanth Nagella, Edoardo Pasca
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.h"
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+*
+* Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+*
+* Output:
+* [1] Regularized image/volume
+
+ * This function is based on the paper by
+* [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+*
+* D. Kazantsev, 2016-18
+*/
+#ifdef __cplusplus
+extern "C" {
+#endif
+//float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float TV_main(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ);
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/compile_mex.m b/Wrappers/Matlab/mex_compile/compile_mex.m
index e1debf3..ee85b49 100644
--- a/Wrappers/Matlab/mex_compile/compile_mex.m
+++ b/Wrappers/Matlab/mex_compile/compile_mex.m
@@ -7,13 +7,14 @@ cd regularizers_CPU/
% compile C regularizers
+mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
mex LLT_model.c LLT_model_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
mex FGP_TV.c FGP_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
mex SplitBregman_TV.c SplitBregman_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
mex TGV_PD.c TGV_PD_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
mex PatchBased_Regul.c PatchBased_Regul_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-delete LLT_model_core.c LLT_model_core.h FGP_TV_core.c FGP_TV_core.h SplitBregman_TV_core.c SplitBregman_TV_core.h TGV_PD_core.c TGV_PD_core.h PatchBased_Regul_core.c PatchBased_Regul_core.h utils.c utils.h CCPiDefines.h
+delete ROF_TV_core.c ROF_TV_core.h LLT_model_core.c LLT_model_core.h FGP_TV_core.c FGP_TV_core.h SplitBregman_TV_core.c SplitBregman_TV_core.h TGV_PD_core.c TGV_PD_core.h PatchBased_Regul_core.c PatchBased_Regul_core.h utils.c utils.h CCPiDefines.h
% compile CUDA-based regularizers
%cd regularizers_GPU/
diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
new file mode 100644
index 0000000..a800add
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
@@ -0,0 +1,106 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "ROF_TV_core.h"
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ * compile: mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+ * D. Kazantsev, 2016-18
+ */
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int i, number_of_dims, iter_numb, dimX, dimY, dimZ;
+ const int *dim_array;
+ float *A, *B, *D1, *D2, *D3, lambda, tau;
+
+ dim_array = mxGetDimensions(prhs[0]);
+ number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+
+ /*Handling Matlab input data*/
+ A = (float *) mxGetData(prhs[0]);
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ tau = (float) mxGetScalar(prhs[2]); /* marching step parameter */
+ iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */
+
+ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+ /*Handling Matlab output data*/
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ /* output arrays*/
+ if (number_of_dims == 2) {
+ dimZ = 0; /*2D case*/
+ B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ D1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ D2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ /* copy into B */
+ copyIm(A, B, dimX, dimY, 1);
+
+ /* start TV iterations */
+ for(i=0; i < iter_numb; i++) {
+
+ /* calculate differences */
+ D1_func(B, D1, dimX, dimY, dimZ);
+ D2_func(B, D2, dimX, dimY, dimZ);
+
+ /* calculate divergence and image update*/
+ TV_main(D1, D2, D2, B, A, lambda, tau, dimX, dimY, dimZ);
+ }
+ }
+
+ if (number_of_dims == 3) {
+ B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ D1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ D2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ D3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+
+ /* copy into B */
+ copyIm(A, B, dimX, dimY, dimZ);
+
+ /* start TV iterations */
+ for(i=0; i < iter_numb; i++) {
+
+ /* calculate differences */
+ D1_func(B, D1, dimX, dimY, dimZ);
+ D2_func(B, D2, dimX, dimY, dimZ);
+ D3_func(B, D3, dimX, dimY, dimZ);
+
+ /* calculate divergence and image update*/
+ TV_main(D1, D2, D3, B, A, lambda, tau, dimX, dimY, dimZ);
+ }
+ }
+
+}
--
cgit v1.2.3