Skip to content

Commit 4fb89fc

Browse files
dkazancdkazanc
authored
Merge pull request #130 from vais-ral/nltv_fix
NLTV fix for nonsquare images
2 parents e2d85c8 + e2601be commit 4fb89fc

File tree

9 files changed

+332
-366
lines changed

9 files changed

+332
-366
lines changed

‎build/run.sh‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,9 @@ rm -r ../build_proj
66
mkdir ../build_proj
77
cd ../build_proj/
88
#make clean
9-
export CIL_VERSION=19.06
9+
export CIL_VERSION=19.09
1010
# install Python modules without CUDA
11-
# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
11+
#cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
1212
# install Python modules with CUDA
1313
cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
1414
# install Matlab modules without CUDA
@@ -21,7 +21,7 @@ make install
2121
cp install/lib/libcilreg.so install/python/ccpi/filters
2222
cd install/python
2323
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
24-
spyder
24+
spyder --new-instance
2525
############### Matlab(linux)###############
2626
# export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6 # if there is libstdc error in matlab
2727
# PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab

‎src/Core/regularisers_CPU/Nonlocal_TV_core.c‎

Lines changed: 39 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
/*
22
* This work is part of the Core Imaging Library developed by
33
* Visual Analytics and Imaging System Group of the Science Technology
4-
* Facilities Council, STFC and Diamond Light Source Ltd.
4+
* Facilities Council, STFC and Diamond Light Source Ltd.
55
*
66
* Copyright 2017 Daniil Kazantsev
77
* Copyright 2017 Srikanth Nagella, Edoardo Pasca
8-
* Copyright 2018 Diamond Light Source Ltd.
8+
* Copyright 2018 Diamond Light Source Ltd.
99
*
1010
* Licensed under the Apache License, Version 2.0 (the "License");
1111
* you may not use this file except in compliance with the License.
@@ -31,60 +31,65 @@
3131
* 2. AR_i - indeces of i neighbours
3232
* 3. AR_j - indeces of j neighbours
3333
* 4. AR_k - indeces of k neighbours (0 - for 2D case)
34-
* 5. Weights_ij(k) - associated weights
34+
* 5. Weights_ij(k) - associated weights
3535
* 6. regularisation parameter
36-
* 7. iterations number
37-
36+
* 7. iterations number
37+
3838
* Output:
39-
* 1. denoised image/volume
39+
* 1. denoised image/volume
4040
* Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060.
41-
41+
4242
*/
4343
/*****************************************************************************/
4444

45-
float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb)
45+
float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb, int switchM)
4646
{
4747

4848
long i, j, k;
4949
int iter;
5050
lambdaReg = 1.0f/lambdaReg;
51-
51+
5252
/*****2D INPUT *****/
5353
if (dimZ == 0) {
5454
copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l);
5555
/* for each pixel store indeces of the most similar neighbours (patches) */
56-
for(iter=0; iter<IterNumb; iter++) {
56+
for(iter=0; iter<IterNumb; iter++) {
5757
#pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, iter) private(i,j)
58+
for(j=0; j<(long)(dimY); j++) {
5859
for(i=0; i<(long)(dimX); i++) {
59-
for(j=0; j<(long)(dimY); j++) {
6060
/*NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);*/ /* NLM - H1 penalty */
61-
NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg); /* NLM - TV penalty */
61+
if (switchM == 1) {
62+
NLM_TV_2D(Output, A_orig, H_j, H_i, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg); /* NLM - TV penalty */
63+
}
64+
else {
65+
NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg); /* NLM - TV penalty */
66+
}
6267
}}
6368
}
64-
}
69+
}
6570
else {
6671
/*****3D INPUT *****/
6772
copyIm(A_orig, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
6873
/* for each pixel store indeces of the most similar neighbours (patches) */
69-
for(iter=0; iter<IterNumb; iter++) {
74+
for(iter=0; iter<IterNumb; iter++) {
7075
#pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, H_k, iter) private(i,j,k)
7176
for(i=0; i<(long)(dimX); i++) {
72-
for(j=0; j<(long)(dimY); j++) {
77+
for(j=0; j<(long)(dimY); j++) {
7378
for(k=0; k<(long)(dimZ); k++) {
7479
/* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambdaReg); */ /* NLM - H1 penalty */
75-
NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambdaReg); /* NLM - TV penalty */
76-
}}}
77-
}
80+
NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambdaReg); /* NLM - TV penalty */
81+
}}}
82+
}
7883
}
7984
return *Output;
8085
}
8186

8287
/***********<<<<Main Function for NLM - H1 penalty>>>>**********/
8388
float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg)
8489
{
85-
long x, i1, j1, index, index_m;
90+
long x, i1, j1, index, index_m;
8691
float value = 0.0f, normweight = 0.0f;
87-
92+
8893
index_m = j*dimX+i;
8994
for(x=0; x < NumNeighb; x++) {
9095
index = (dimX*dimY*x) + j*dimX+i;
@@ -99,17 +104,17 @@ float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_
99104
/*3D version*/
100105
float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg)
101106
{
102-
long x, i1, j1, k1, index;
107+
long x, i1, j1, k1, index;
103108
float value = 0.0f, normweight = 0.0f;
104-
109+
105110
for(x=0; x < NumNeighb; x++) {
106111
index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
107112
i1 = H_i[index];
108113
j1 = H_j[index];
109114
k1 = H_k[index];
110115
value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index];
111116
normweight += Weights[index];
112-
}
117+
}
113118
A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight);
114119
return *A;
115120
}
@@ -118,56 +123,56 @@ float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_
118123
/***********<<<<Main Function for NLM - TV penalty>>>>**********/
119124
float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg)
120125
{
121-
long x, i1, j1, index, index_m;
126+
long x, i1, j1, index, index_m;
122127
float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
123-
128+
124129
index_m = j*dimX+i;
125-
130+
126131
for(x=0; x < NumNeighb; x++) {
127132
index = (dimX*dimY*x) + j*dimX+i; /*c*/
128133
i1 = H_i[index];
129134
j1 = H_j[index];
130135
NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index];
131136
}
132-
137+
133138
NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */
134139
NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS));
135-
140+
136141
for(x=0; x < NumNeighb; x++) {
137142
index = (dimX*dimY*x) + j*dimX+i; /*c*/
138143
i1 = H_i[index];
139144
j1 = H_j[index];
140145
value += A[j1*dimX+i1]*NLCoeff*Weights[index];
141146
normweight += Weights[index]*NLCoeff;
142-
}
147+
}
143148
A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight);
144149
return *A;
145150
}
146151
/*3D version*/
147152
float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg)
148153
{
149-
long x, i1, j1, k1, index;
154+
long x, i1, j1, k1, index;
150155
float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
151-
156+
152157
for(x=0; x < NumNeighb; x++) {
153158
index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
154159
i1 = H_i[index];
155160
j1 = H_j[index];
156161
k1 = H_k[index];
157162
NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index];
158163
}
159-
164+
160165
NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */
161166
NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS));
162-
167+
163168
for(x=0; x < NumNeighb; x++) {
164169
index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
165170
i1 = H_i[index];
166171
j1 = H_j[index];
167172
k1 = H_k[index];
168173
value += A[(dimX*dimY*k1) + j1*dimX+i1]*NLCoeff*Weights[index];
169174
normweight += Weights[index]*NLCoeff;
170-
}
175+
}
171176
A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight);
172177
return *A;
173178
}

‎src/Core/regularisers_CPU/Nonlocal_TV_core.h‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
/*
22
* This work is part of the Core Imaging Library developed by
33
* Visual Analytics and Imaging System Group of the Science Technology
4-
* Facilities Council, STFC and Diamond Light Source Ltd.
4+
* Facilities Council, STFC and Diamond Light Source Ltd.
55
*
66
* Copyright 2017 Daniil Kazantsev
77
* Copyright 2017 Srikanth Nagella, Edoardo Pasca
8-
* Copyright 2018 Diamond Light Source Ltd.
8+
* Copyright 2018 Diamond Light Source Ltd.
99
*
1010
* Licensed under the Apache License, Version 2.0 (the "License");
1111
* you may not use this file except in compliance with the License.
@@ -39,19 +39,19 @@
3939
* 2. AR_i - indeces of i neighbours
4040
* 3. AR_j - indeces of j neighbours
4141
* 4. AR_k - indeces of k neighbours (0 - for 2D case)
42-
* 5. Weights_ij(k) - associated weights
42+
* 5. Weights_ij(k) - associated weights
4343
* 6. regularisation parameter
44-
* 7. iterations number
45-
44+
* 7. iterations number
45+
4646
* Output:
47-
* 1. denoised image/volume
48-
* Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060.
47+
* 1. denoised image/volume
48+
* Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060.
4949
*/
50-
50+
5151
#ifdef __cplusplus
5252
extern "C" {
5353
#endif
54-
CCPI_EXPORT float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb);
54+
CCPI_EXPORT float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb, int switchM);
5555
CCPI_EXPORT float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg);
5656
CCPI_EXPORT float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg);
5757
CCPI_EXPORT float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg);

0 commit comments

Comments
 (0)