summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Core/CMakeLists.txt31
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.c76
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.h11
-rw-r--r--src/Core/regularisers_CPU/PD_TV_core.c166
-rw-r--r--src/Core/regularisers_CPU/PD_TV_core.h57
-rw-r--r--src/Core/regularisers_CPU/utils.c46
-rw-r--r--src/Core/regularisers_CPU/utils.h1
-rw-r--r--src/Matlab/mex_compile/compileCPU_mex_Linux.m7
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/PD_TV.c98
-rw-r--r--src/Python/ccpi/filters/regularisers.py28
-rw-r--r--src/Python/setup-regularisers.py.in25
-rw-r--r--src/Python/src/cpu_regularisers.pyx42
12 files changed, 489 insertions, 99 deletions
diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index eea0d63..ac7ec91 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -20,7 +20,7 @@ if (OPENMP_FOUND)
set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
-
+
endif()
## Build the regularisers package as a library
@@ -39,21 +39,21 @@ if(WIN32)
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib")
-
+
set (EXTRA_LIBRARIES)
-
+
message("library lib: ${LIBRARY_LIB}")
-
+
elseif(UNIX)
- set (FLAGS "-O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS ")
+ set (FLAGS "-O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS ")
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
-
- set (EXTRA_LIBRARIES
+
+ set (EXTRA_LIBRARIES
"gomp"
"m"
)
-
+
endif()
message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
@@ -66,6 +66,7 @@ message("Adding regularisers as a shared library")
add_library(cilreg SHARED
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c
@@ -80,8 +81,8 @@ add_library(cilreg SHARED
${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/NonlocalMarching_Inpaint_core.c
)
target_link_libraries(cilreg ${EXTRA_LIBRARIES} )
-include_directories(cilreg PUBLIC
- ${LIBRARY_INC}/include
+include_directories(cilreg PUBLIC
+ ${LIBRARY_INC}/include
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/
${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/ )
@@ -92,14 +93,14 @@ if (UNIX)
message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
install(TARGETS cilreg
LIBRARY DESTINATION lib
- CONFIGURATIONS ${CMAKE_BUILD_TYPE}
+ CONFIGURATIONS ${CMAKE_BUILD_TYPE}
)
elseif(WIN32)
message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin")
- install(TARGETS cilreg
+ install(TARGETS cilreg
RUNTIME DESTINATION bin
ARCHIVE DESTINATION lib
- CONFIGURATIONS ${CMAKE_BUILD_TYPE}
+ CONFIGURATIONS ${CMAKE_BUILD_TYPE}
)
endif()
@@ -126,14 +127,14 @@ if (BUILD_CUDA)
message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
install(TARGETS cilregcuda
LIBRARY DESTINATION lib
- CONFIGURATIONS ${CMAKE_BUILD_TYPE}
+ CONFIGURATIONS ${CMAKE_BUILD_TYPE}
)
elseif(WIN32)
message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin")
install(TARGETS cilregcuda
RUNTIME DESTINATION bin
ARCHIVE DESTINATION lib
- CONFIGURATIONS ${CMAKE_BUILD_TYPE}
+ CONFIGURATIONS ${CMAKE_BUILD_TYPE}
)
endif()
else()
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index a16a2e5..ff67af2 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -46,12 +46,12 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
float tk = 1.0f;
float tkp1 =1.0f;
int count = 0;
-
+
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);
-
+
if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
P2 = calloc(DimTotal, sizeof(float));
@@ -59,32 +59,32 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
P2_prev = calloc(DimTotal, sizeof(float));
R1 = calloc(DimTotal, sizeof(float));
R2 = calloc(DimTotal, sizeof(float));
-
+
/* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
-
+
if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
/* computing the gradient of the objective function */
Obj_func2D(Input, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+
/* apply nonnegativity */
if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+
/*Taking a step towards minus of the gradient*/
Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+
/* projection step */
Proj_func2D(P1, P2, methodTV, DimTotal);
-
+
/*updating R and t*/
tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-
+
/*storing old values*/
copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
tk = tkp1;
-
+
/* check early stopping criteria */
if ((epsil != 0.0f) && (ll % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
@@ -105,7 +105,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
/*3D case*/
float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=NULL;
DimTotal = (long)(dimX*dimY*dimZ);
-
+
if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
P2 = calloc(DimTotal, sizeof(float));
@@ -116,28 +116,28 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
R1 = calloc(DimTotal, sizeof(float));
R2 = calloc(DimTotal, sizeof(float));
R3 = calloc(DimTotal, sizeof(float));
-
+
/* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
-
+
if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* computing the gradient of the objective function */
Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* apply nonnegativity */
if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+
/*Taking a step towards minus of the gradient*/
Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* projection step */
Proj_func3D(P1, P2, P3, methodTV, DimTotal);
-
+
/*updating R and t*/
tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-
+
/* calculate norm - stopping rules*/
if ((epsil != 0.0f) && (ll % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
@@ -151,22 +151,22 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
if (re < epsil) count++;
if (count > 3) break;
}
-
+
/*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 (epsil != 0.0f) 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 */
-
+
return 0;
}
@@ -202,36 +202,6 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la
}}
return 1;
}
-float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
-{
- float val1, val2, denom, sq_denom;
- long i;
- if (methTV == 0) {
- /* isotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
- for(i=0; i<DimTotal; i++) {
- denom = powf(P1[i],2) + powf(P2[i],2);
- if (denom > 1.0f) {
- sq_denom = 1.0f/sqrtf(denom);
- P1[i] = P1[i]*sq_denom;
- P2[i] = P2[i]*sq_denom;
- }
- }
- }
- else {
- /* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(i,val1,val2)
- for(i=0; i<DimTotal; i++) {
- val1 = fabs(P1[i]);
- val2 = fabs(P2[i]);
- if (val1 < 1.0f) {val1 = 1.0f;}
- if (val2 < 1.0f) {val2 = 1.0f;}
- P1[i] = P1[i]/val1;
- P2[i] = P2[i]/val2;
- }
- }
- return 1;
-}
float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal)
{
long i;
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h
index 04e6e80..4466a35 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -29,12 +29,12 @@ limitations under the License.
/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
*
* Input Parameters:
- * 1. Noisy image/volume
- * 2. lambda - regularization parameter
+ * 1. Noisy image/volume
+ * 2. lambda - regularization parameter
* 3. Number of iterations
- * 4. eplsilon: tolerance constant
+ * 4. eplsilon: tolerance constant
* 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
- * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
*
* Output:
* [1] Filtered/regularized image/volume
@@ -44,7 +44,7 @@ limitations under the License.
* 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"
*/
-
+
#ifdef __cplusplus
extern "C" {
#endif
@@ -52,7 +52,6 @@ CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float *infovector
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);
-CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal);
CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal);
CCPI_EXPORT float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c
new file mode 100644
index 0000000..a8c3cfb
--- /dev/null
+++ b/src/Core/regularisers_CPU/PD_TV_core.c
@@ -0,0 +1,166 @@
+/*
+ * 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 2019 Daniil Kazantsev
+ * Copyright 2019 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 "PD_TV_core.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. lipschitz_const: convergence related parameter
+ * 6. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 7. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
+
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
+ */
+
+float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
+{
+ int ll;
+ long j, DimTotal;
+ float re, re1, tau, sigma, theta, lt;
+ re = 0.0f; re1 = 0.0f;
+ int count = 0;
+
+ tau = 1.0/powf(lipschitz_const,0.5);
+ sigma = 1.0/powf(lipschitz_const,0.5);
+ theta = 1.0f;
+ lt = tau/lambdaPar;
+ ll = 0;
+
+ copyIm(Input, U, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+ if (dimZ <= 1) {
+ /*2D case */
+ float *U_old=NULL, *P1=NULL, *P2=NULL;
+ DimTotal = (long)(dimX*dimY);
+
+ U_old = calloc(DimTotal, sizeof(float));
+ P1 = calloc(DimTotal, sizeof(float));
+ P2 = calloc(DimTotal, sizeof(float));
+
+ /* begin iterations */
+ for(ll=0; ll<iterationsNumb; ll++) {
+
+ //if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
+ /* computing the gradient of the objective function */
+ DualP2D(U, P1, P2, (long)(dimX), (long)(dimY), sigma);
+
+ /* apply nonnegativity */
+ if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (U[j] < 0.0f) U[j] = 0.0f;}
+
+ /* projection step */
+ Proj_func2D(P1, P2, methodTV, DimTotal);
+
+ /* copy U to U_old */
+ copyIm(U, U_old, (long)(dimX), (long)(dimY), 1l);
+
+ DivProj2D(U, Input, P1, P2,(long)(dimX), (long)(dimY), lt, tau);
+
+ /* check early stopping criteria */
+ if ((epsil != 0.0f) && (ll % 5 == 0)) {
+ re = 0.0f; re1 = 0.0f;
+ for(j=0; j<DimTotal; j++)
+ {
+ re += powf(U[j] - U_old[j],2);
+ re1 += powf(U[j],2);
+ }
+ re = sqrtf(re)/sqrtf(re1);
+ if (re < epsil) count++;
+ if (count > 3) break;
+ }
+ /*get updated solution*/
+
+ getX2D(U, U_old, theta, DimTotal);
+ }
+ free(P1); free(P2); free(U_old);
+ }
+ else {
+ /*3D case*/
+ }
+ /*adding info into info_vector */
+ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
+ infovector[1] = re; /* reached tolerance */
+
+ return 0;
+}
+
+/*****************************************************************/
+/************************2D-case related Functions */
+/*****************************************************************/
+
+/*Calculating dual variable (using forward differences)*/
+float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma)
+{
+ long i,j,index;
+ #pragma omp parallel for shared(U,P1,P2) private(index,i,j)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ index = j*dimX+i;
+ /* symmetric boundary conditions (Neuman) */
+ if (i == dimX-1) P1[index] += sigma*(U[j*dimX+(i-1)] - U[index]);
+ else P1[index] += sigma*(U[j*dimX+(i+1)] - U[index]);
+ if (j == dimY-1) P2[index] += sigma*(U[(j-1)*dimX+i] - U[index]);
+ else P2[index] += sigma*(U[(j+1)*dimX+i] - U[index]);
+ }}
+ return 1;
+}
+
+/* Divergence for P dual */
+float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau)
+{
+ long i,j,index;
+ float P_v1, P_v2, div_var;
+ #pragma omp parallel for shared(U,Input,P1,P2) private(i, j, P_v1, P_v2, div_var)
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ index = j*dimX+i;
+ /* symmetric boundary conditions (Neuman) */
+ if (i == 0) P_v1 = -P1[index];
+ else P_v1 = -(P1[index] - P1[j*dimX+(i-1)]);
+ if (j == 0) P_v2 = -P2[index];
+ else P_v2 = -(P2[index] - P2[(j-1)*dimX+i]);
+ div_var = P_v1 + P_v2;
+ U[index] = (U[index] - tau*div_var + lt*Input[index])/(1.0 + lt);
+ }}
+ return *U;
+}
+
+/*get the updated solution*/
+float getX2D(float *U, float *U_old, float theta, long DimTotal)
+{
+ long i;
+ #pragma omp parallel for shared(U,U_old) private(i)
+ for(i=0; i<DimTotal; i++) {
+ U[i] += theta*(U[i] - U_old[i]);
+ }
+ return *U;
+}
+
+
+/*****************************************************************/
+/************************3D-case related Functions */
+/*****************************************************************/
diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h
new file mode 100644
index 0000000..b52689a
--- /dev/null
+++ b/src/Core/regularisers_CPU/PD_TV_core.h
@@ -0,0 +1,57 @@
+/*
+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 2019 Daniil Kazantsev
+Copyright 2019 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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. lipschitz_const: convergence related parameter
+ * 6. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 7. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
+
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
+
+CCPI_EXPORT float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma);
+CCPI_EXPORT float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau);
+CCPI_EXPORT float getX2D(float *U, float *U_old, float theta, long DimTotal);
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c
index 5bb3a5c..cf4ad72 100644
--- a/src/Core/regularisers_CPU/utils.c
+++ b/src/Core/regularisers_CPU/utils.c
@@ -65,7 +65,7 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int
{
int i, j, i1, j1, index;
float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
-
+
/* first calculate \grad U_xy*/
for(j=0; j<dimY; j++) {
for(i=0; i<dimX; i++) {
@@ -73,7 +73,7 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int
/* boundary conditions */
i1 = i + 1; if (i == dimX-1) i1 = i;
j1 = j + 1; if (j == dimY-1) j1 = j;
-
+
/* 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+ */
@@ -90,7 +90,7 @@ float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int
{
long i, j, k, i1, j1, k1, index;
float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f;
-
+
/* first calculate \grad U_xy*/
for(j=0; j<(long)(dimY); j++) {
for(i=0; i<(long)(dimX); i++) {
@@ -100,12 +100,12 @@ float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int
i1 = i + 1; if (i == (long)(dimX-1)) i1 = i;
j1 = j + 1; if (j == (long)(dimY-1)) j1 = j;
k1 = k + 1; if (k == (long)(dimZ-1)) k1 = k;
-
+
/* Forward differences */
NOMx_2 = powf((float)(U[(dimX*dimY)*k + j1*dimX+i] - U[index]),2); /* x+ */
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 += 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 */
}
@@ -131,12 +131,12 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
x_diff = (x_ratio * j) - x;
y_diff = (y_ratio * i) - y;
index = y*w+x ;
-
+
A = Input[index];
B = Input[index+1];
C = Input[index+w];
D = Input[index+w+1];
-
+
gray = (float)(A*(1.0 - x_diff)*(1.0 - y_diff) + B*(x_diff)*(1.0 - y_diff) +
C*(y_diff)*(1.0 - x_diff) + D*(x_diff*y_diff));
@@ -144,3 +144,35 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
}}
return *Scaled;
}
+
+/*2D Projection onto convex set for P*/
+float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
+{
+ float val1, val2, denom, sq_denom;
+ long i;
+ if (methTV == 0) {
+ /* isotropic TV*/
+#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
+ for(i=0; i<DimTotal; i++) {
+ denom = powf(P1[i],2) + powf(P2[i],2);
+ if (denom > 1.0f) {
+ sq_denom = 1.0f/sqrtf(denom);
+ P1[i] = P1[i]*sq_denom;
+ P2[i] = P2[i]*sq_denom;
+ }
+ }
+ }
+ else {
+ /* anisotropic TV*/
+#pragma omp parallel for shared(P1,P2) private(i,val1,val2)
+ for(i=0; i<DimTotal; i++) {
+ val1 = fabs(P1[i]);
+ val2 = fabs(P2[i]);
+ if (val1 < 1.0f) {val1 = 1.0f;}
+ if (val2 < 1.0f) {val2 = 1.0f;}
+ P1[i] = P1[i]/val1;
+ P2[i] = P2[i]/val2;
+ }
+ }
+ return 1;
+}
diff --git a/src/Core/regularisers_CPU/utils.h b/src/Core/regularisers_CPU/utils.h
index 8f0ba59..e050f86 100644
--- a/src/Core/regularisers_CPU/utils.h
+++ b/src/Core/regularisers_CPU/utils.h
@@ -31,6 +31,7 @@ CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, i
CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
CCPI_EXPORT float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2);
+CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal);
#ifdef __cplusplus
}
#endif
diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 2ed7ea6..5330d7f 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -28,6 +28,10 @@ movefile('FGP_TV.mex*',Pathmove);
fprintf('%s \n', 'Compiling SB-TV...');
mex SB_TV.c SB_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
movefile('SB_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling PD-TV...');
+mex PD_TV.c PD_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('PD_TV.mex*',Pathmove);
fprintf('%s \n', 'Compiling dFGP-TV...');
mex FGP_dTV.c FGP_dTV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
@@ -75,7 +79,8 @@ movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* LLT_ROF_core* CCPiDefines.h
delete PatchSelect_core* Nonlocal_TV_core*
delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
-fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
+delete PD_TV_core*
+fprintf('%s \n', '<<<<<<< CPU regularisers were successfully compiled! >>>>>>>');
pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i);
cd(pathA2); \ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
new file mode 100644
index 0000000..eac2d18
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
@@ -0,0 +1,98 @@
+/*
+ * 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 2019 Daniil Kazantsev
+ * Copyright 2019 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 "PD_TV_core.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
+ * 7. lipschitz_const: convergence related parameter
+
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
+ */
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+
+{
+ int number_of_dims, iter, methTV, nonneg;
+ mwSize dimX, dimY, dimZ;
+ const mwSize *dim_array;
+ float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const;
+
+ 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, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const");
+
+ Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
+ iter = 400; /* default iterations number */
+ epsil = 1.0e-06; /* default tolerance constant */
+ methTV = 0; /* default isotropic TV penalty */
+ nonneg = 0; /* default nonnegativity switch, off - 0 */
+ lipschitz_const = 12.0; /* lipschitz_const */
+
+ 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)) {
+ 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)) {
+ 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) lipschitz_const = (float) mxGetScalar(prhs[6]);
+
+
+ /*Handling Matlab output data*/
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ if (number_of_dims == 2) {
+ dimZ = 1; /*2D case*/
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ if (number_of_dims == 3) {
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ }
+ mwSize vecdim[1];
+ vecdim[0] = 2;
+ infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));
+
+ /* running the function */
+ PDTV_CPU_main(Input, Output, infovec, lambda, iter, epsil, lipschitz_const, methTV, nonneg, dimX, dimY, dimZ);
+}
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 0b5b2ee..d65c0b9 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -2,7 +2,7 @@
script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
"""
-from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_PD_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
try:
from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU
gpu_enabled = True
@@ -51,6 +51,31 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
raise ValueError ('GPU is not available')
raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
.format(device))
+
+def PD_TV(inputData, regularisation_parameter, iterations,
+ tolerance_param, methodTV, nonneg, lipschitz_const, device='cpu'):
+ if device == 'cpu':
+ return TV_PD_CPU(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ lipschitz_const)
+ elif device == 'gpu' and gpu_enabled:
+ return TV_PD_CPU(inputData,
+ regularisation_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ lipschitz_const)
+ else:
+ if not gpu_enabled and device == 'gpu':
+ raise ValueError ('GPU is not available')
+ raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+ .format(device))
+
def SB_TV(inputData, regularisation_parameter, iterations,
tolerance_param, methodTV, device='cpu'):
if device == 'cpu':
@@ -212,4 +237,3 @@ def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, itera
def NVM_INP(inputData, maskData, SW_increment, iterations):
return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
-
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 4c578e3..9bcd46d 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -8,13 +8,13 @@ from Cython.Distutils import build_ext
import os
import sys
import numpy
-import platform
+import platform
cil_version=os.environ['CIL_VERSION']
if cil_version == '':
print("Please set the environmental variable CIL_VERSION")
sys.exit(1)
-
+
library_include_path = ""
library_lib_path = ""
try:
@@ -23,7 +23,7 @@ try:
except:
library_include_path = os.environ['PREFIX']+'/include'
pass
-
+
extra_include_dirs = [numpy.get_include(), library_include_path]
#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
extra_compile_args = []
@@ -38,6 +38,7 @@ extra_include_dirs += [os.path.join(".." , "Core"),
os.path.join(".." , "Core", "regularisers_CPU"),
os.path.join(".." , "Core", "inpainters_CPU"),
os.path.join(".." , "Core", "regularisers_GPU" , "TV_FGP" ) ,
+ os.path.join(".." , "Core", "regularisers_GPU" , "TV_PD" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TV_ROF" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TV_SB" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TGV" ) ,
@@ -48,12 +49,12 @@ extra_include_dirs += [os.path.join(".." , "Core"),
os.path.join(".." , "Core", "regularisers_GPU" , "PatchSelect" ) ,
"."]
-if platform.system() == 'Windows':
- extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]
+if platform.system() == 'Windows':
+ extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]
else:
extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
extra_libraries += [@EXTRA_OMP_LIB@]
-
+
setup(
name='ccpi',
description='CCPi Core Imaging Library - Image regularisers',
@@ -61,13 +62,13 @@ setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("ccpi.filters.cpu_regularisers",
sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
- include_dirs=extra_include_dirs,
- library_dirs=extra_library_dirs,
- extra_compile_args=extra_compile_args,
- libraries=extra_libraries ),
-
+ include_dirs=extra_include_dirs,
+ library_dirs=extra_library_dirs,
+ extra_compile_args=extra_compile_args,
+ libraries=extra_libraries ),
+
],
- zip_safe = False,
+ zip_safe = False,
packages = {'ccpi', 'ccpi.filters', 'ccpi.supp'},
)
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 4917d06..724634b 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -20,6 +20,7 @@ cimport numpy as np
cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iterationsNumb, float tau, float epsil, 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 PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ);
cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
@@ -58,9 +59,6 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
np.ones([2], dtype='float32')
- # Run ROF iterations for 2D data
- # TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)
- # Run ROF iterations for 2D data
if isinstance (regularisation_parameter, np.ndarray):
reg = regularisation_parameter.copy()
TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &reg[0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)
@@ -158,6 +156,44 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
dims[2], dims[1], dims[0])
return (outputData,infovec)
+#****************************************************************#
+#****************** Total-variation Primal-dual *****************#
+#****************************************************************#
+def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const):
+ if inputData.ndim == 2:
+ return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
+ elif inputData.ndim == 3:
+ return 0
+
+def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterationsNumb,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ float lipschitz_const):
+
+ 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=1, mode="c"] infovec = \
+ np.ones([2], dtype='float32')
+
+ #/* Run FGP-TV iterations for 2D data */
+ PDTV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter,
+ iterationsNumb,
+ tolerance_param,
+ lipschitz_const,
+ methodTV,
+ nonneg,
+ dims[1],dims[0], 1)
+
+ return (outputData,infovec)
+
#***************************************************************#
#********************** Total-variation SB *********************#
#***************************************************************#
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-09 23:02:43 +0000