Skip to content

Commit

Permalink
Merge pull request #128 from vais-ral/fgp_cpu_resort_for
Browse files Browse the repository at this point in the history 
Fgp cpu resort for
  • Loading branch information
@dkazanc
dkazanc authored Sep 26, 2019
2 parents 4fb89fc + 4c67694 commit 51e49a4
Show file tree
Hide file tree
Showing 12 changed files with 1,288 additions and 1,295 deletions.
29 changes: 14 additions & 15 deletions demos/Matlab_demos/demoMatlab_3Ddenoise.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,19 +182,18 @@
tic; u_fgp_dtv = FGP_dTV(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc;
figure; imshow(u_fgp_dtv(:,:,7), [0 1]); title('FGP-dTV denoised volume (CPU)');
%%
fprintf('Denoise a volume using the FGP-dTV model (GPU) \n');

% create another volume (reference) with slightly less amount of noise
vol3D_ref = zeros(N,N,slices, 'single');
for i = 1:slices
vol3D_ref(:,:,i) = Im + .01*randn(size(Im));
end
vol3D_ref(vol3D_ref < 0) = 0;
% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV)

iter_fgp = 300; % number of FGP iterations
epsil_tol = 0.0; % tolerance
eta = 0.2; % Reference image gradient smoothing constant
tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc;
figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)');
% fprintf('Denoise a volume using the FGP-dTV model (GPU) \n');
% % create another volume (reference) with slightly less amount of noise
% vol3D_ref = zeros(N,N,slices, 'single');
% for i = 1:slices
% vol3D_ref(:,:,i) = Im + .01*randn(size(Im));
% end
% vol3D_ref(vol3D_ref < 0) = 0;
% % vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV)
%
% iter_fgp = 300; % number of FGP iterations
% epsil_tol = 0.0; % tolerance
% eta = 0.2; % Reference image gradient smoothing constant
% tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc;
% figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)');
%%
1 change: 0 additions & 1 deletion demos/Matlab_demos/demoMatlab_denoise.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,6 @@
figure; imagesc(u_nltv, [0 1]); colormap(gray); daspect([1 1 1]); title('Non-local Total Variation denoised image (CPU)');
%%
%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< %

fprintf('Denoise using the FGP-dTV model (CPU) \n');
% create another image (reference) with slightly less amount of noise
u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0;
Expand Down
283 changes: 137 additions & 146 deletions src/Core/regularisers_CPU/Diffus4th_order_core.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,35 +50,35 @@ float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float l
float *W_Lapl=NULL, *Output_prev=NULL;
sigmaPar2 = sigmaPar*sigmaPar;
DimTotal = dimX*dimY*dimZ;

W_Lapl = calloc(DimTotal, sizeof(float));

if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));

/* copy into output */
copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));

for(i=0; i < iterationsNumb; i++) {
if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));

if (dimZ == 1) {
if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
if (dimZ == 1) {
/* running 2D diffusion iterations */
/* Calculating weighted Laplacian */
Weighted_Laplc2D(W_Lapl, Output, sigmaPar2, dimX, dimY);
/* Perform iteration step */
Diffusion_update_step2D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY));
}
else {
}
else {
/* running 3D diffusion iterations */
/* Calculating weighted Laplacian */
Weighted_Laplc3D(W_Lapl, Output, sigmaPar2, dimX, dimY, dimZ);
/* Perform iteration step */
Diffusion_update_step3D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
}

/* check early stopping criteria */
if ((epsil != 0.0f) && (i % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
}
/* check early stopping criteria */
if ((epsil != 0.0f) && (i % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
for(j=0; j<DimTotal; j++)
{
re += powf(Output[j] - Output_prev[j],2);
Expand All @@ -87,10 +87,10 @@ float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float l
re = sqrtf(re)/sqrtf(re1);
if (re < epsil) count++;
if (count > 3) break;
}
}
free(W_Lapl);

}
}
free(W_Lapl);
if (epsil != 0.0f) free(Output_prev);
/*adding info into info_vector */
infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/
Expand All @@ -104,74 +104,71 @@ float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, long dimX, long di
{
long i,j,i1,i2,j1,j2,index;
float gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq;

#pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq)

#pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq)
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;

index = j*dimX+i;

gradX = 0.5f*(U0[j*dimX+i2] - U0[j*dimX+i1]);
gradX_sq = pow(gradX,2);

gradY = 0.5f*(U0[j2*dimX+i] - U0[j1*dimX+i]);
gradY_sq = pow(gradY,2);

gradXX = U0[j*dimX+i2] + U0[j*dimX+i1] - 2*U0[index];
gradYY = U0[j2*dimX+i] + U0[j1*dimX+i] - 2*U0[index];

gradXY = 0.25f*(U0[j2*dimX+i2] + U0[j1*dimX+i1] - U0[j1*dimX+i2] - U0[j2*dimX+i1]);
xy_2 = 2.0f*gradX*gradY*gradXY;

denom = gradX_sq + gradY_sq;

if (denom <= EPS) {
V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/EPS;
V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/EPS;
}
else {
V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/denom;
V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/denom;
}

c = 1.0f/(1.0f + denom/sigma);
c_sq = c*c;

W_Lapl[index] = c_sq*V_norm + c*V_orth;
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
index = j*dimX+i;

gradX = 0.5f*(U0[j*dimX+i2] - U0[j*dimX+i1]);
gradX_sq = pow(gradX,2);

gradY = 0.5f*(U0[j2*dimX+i] - U0[j1*dimX+i]);
gradY_sq = pow(gradY,2);

gradXX = U0[j*dimX+i2] + U0[j*dimX+i1] - 2*U0[index];
gradYY = U0[j2*dimX+i] + U0[j1*dimX+i] - 2*U0[index];

gradXY = 0.25f*(U0[j2*dimX+i2] + U0[j1*dimX+i1] - U0[j1*dimX+i2] - U0[j2*dimX+i1]);
xy_2 = 2.0f*gradX*gradY*gradXY;

denom = gradX_sq + gradY_sq;

if (denom <= EPS) {
V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/EPS;
V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/EPS;
}
}
return *W_Lapl;
else {
V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/denom;
V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/denom;
}

c = 1.0f/(1.0f + denom/sigma);
c_sq = c*c;

W_Lapl[index] = c_sq*V_norm + c*V_orth;
}}
return *W_Lapl;
}

float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY)
{
long i,j,i1,i2,j1,j2,index;
long i,j,i1,i2,j1,j2,index;
float gradXXc, gradYYc;

#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,index,gradXXc,gradYYc)

#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,index,gradXXc,gradYYc)
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;
index = j*dimX+i;

gradXXc = W_Lapl[j*dimX+i2] + W_Lapl[j*dimX+i1] - 2*W_Lapl[index];
gradYYc = W_Lapl[j2*dimX+i] + W_Lapl[j1*dimX+i] - 2*W_Lapl[index];

Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc) - (Output[index] - Input[index]));
}
}
return *Output;
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
index = j*dimX+i;

gradXXc = W_Lapl[j*dimX+i2] + W_Lapl[j*dimX+i1] - 2*W_Lapl[index];
gradYYc = W_Lapl[j2*dimX+i] + W_Lapl[j1*dimX+i] - 2*W_Lapl[index];

Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc) - (Output[index] - Input[index]));
}}
return *Output;
}
/********************************************************************/
/***************************3D Functions*****************************/
Expand All @@ -180,95 +177,89 @@ float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, long dimX, long di
{
long i,j,k,i1,i2,j1,j2,k1,k2,index;
float gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq, gradZ, gradZ_sq, gradZZ, gradXZ, gradYZ, xyz_1, xyz_2;

#pragma omp parallel for shared(W_Lapl) private(i,j,k,i1,i2,j1,j2,k1,k2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq, gradZ, gradZ_sq, gradZZ, gradXZ, gradYZ, xyz_1, xyz_2)
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;

for(k=0; k<dimZ; k++) {
/* symmetric boundary conditions */
k1 = k+1; if (k1 == dimZ) k1 = k-1;
k2 = k-1; if (k2 < 0) k2 = k+1;

index = (dimX*dimY)*k + j*dimX+i;

gradX = 0.5f*(U0[(dimX*dimY)*k + j*dimX+i2] - U0[(dimX*dimY)*k + j*dimX+i1]);
gradX_sq = pow(gradX,2);

gradY = 0.5f*(U0[(dimX*dimY)*k + j2*dimX+i] - U0[(dimX*dimY)*k + j1*dimX+i]);

#pragma omp parallel for shared(W_Lapl) private(i,j,k,i1,i2,j1,j2,k1,k2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq, gradZ, gradZ_sq, gradZZ, gradXZ, gradYZ, xyz_1, xyz_2)
for(k=0; k<dimZ; k++) {
/* symmetric boundary conditions */
k1 = k+1; if (k1 == dimZ) k1 = k-1;
k2 = k-1; if (k2 < 0) k2 = k+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;

index = (dimX*dimY)*k + j*dimX+i;

gradX = 0.5f*(U0[(dimX*dimY)*k + j*dimX+i2] - U0[(dimX*dimY)*k + j*dimX+i1]);
gradX_sq = pow(gradX,2);

gradY = 0.5f*(U0[(dimX*dimY)*k + j2*dimX+i] - U0[(dimX*dimY)*k + j1*dimX+i]);
gradY_sq = pow(gradY,2);

gradZ = 0.5f*(U0[(dimX*dimY)*k2 + j*dimX+i] - U0[(dimX*dimY)*k1 + j*dimX+i]);
gradZ_sq = pow(gradZ,2);

gradXX = U0[(dimX*dimY)*k + j*dimX+i2] + U0[(dimX*dimY)*k + j*dimX+i1] - 2*U0[index];
gradYY = U0[(dimX*dimY)*k + j2*dimX+i] + U0[(dimX*dimY)*k + j1*dimX+i] - 2*U0[index];
gradZZ = U0[(dimX*dimY)*k2 + j*dimX+i] + U0[(dimX*dimY)*k1 + j*dimX+i] - 2*U0[index];

gradXY = 0.25f*(U0[(dimX*dimY)*k + j2*dimX+i2] + U0[(dimX*dimY)*k + j1*dimX+i1] - U0[(dimX*dimY)*k + j1*dimX+i2] - U0[(dimX*dimY)*k + j2*dimX+i1]);
gradXZ = 0.25f*(U0[(dimX*dimY)*k2 + j*dimX+i2] - U0[(dimX*dimY)*k2+j*dimX+i1] - U0[(dimX*dimY)*k1+j*dimX+i2] + U0[(dimX*dimY)*k1+j*dimX+i1]);
gradYZ = 0.25f*(U0[(dimX*dimY)*k2 +j2*dimX+i] - U0[(dimX*dimY)*k2+j1*dimX+i] - U0[(dimX*dimY)*k1+j2*dimX+i] + U0[(dimX*dimY)*k1+j1*dimX+i]);

xy_2 = 2.0f*gradX*gradY*gradXY;
xyz_1 = 2.0f*gradX*gradZ*gradXZ;
xyz_2 = 2.0f*gradY*gradZ*gradYZ;

denom = gradX_sq + gradY_sq + gradZ_sq;

if (denom <= EPS) {
V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/EPS;
if (denom <= EPS) {
V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/EPS;
V_orth = ((gradY_sq + gradZ_sq)*gradXX + (gradX_sq + gradZ_sq)*gradYY + (gradX_sq + gradY_sq)*gradZZ - xy_2 - xyz_1 - xyz_2)/EPS;
}
else {
V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/denom;
}
else {
V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/denom;
V_orth = ((gradY_sq + gradZ_sq)*gradXX + (gradX_sq + gradZ_sq)*gradYY + (gradX_sq + gradY_sq)*gradZZ - xy_2 - xyz_1 - xyz_2)/denom;
}

}
c = 1.0f/(1.0f + denom/sigma);
c_sq = c*c;

W_Lapl[index] = c_sq*V_norm + c*V_orth;
}
}
}
return *W_Lapl;
}}}
return *W_Lapl;
}

float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY, long dimZ)
{
long i,j,i1,i2,j1,j2,index,k,k1,k2;
long i,j,i1,i2,j1,j2,index,k,k1,k2;
float gradXXc, gradYYc, gradZZc;

#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,k,k1,k2,index,gradXXc,gradYYc,gradZZc)
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;

for(k=0; k<dimZ; k++) {
/* symmetric boundary conditions */
k1 = k+1; if (k1 == dimZ) k1 = k-1;
k2 = k-1; if (k2 < 0) k2 = k+1;

index = (dimX*dimY)*k + j*dimX+i;

gradXXc = W_Lapl[(dimX*dimY)*k + j*dimX+i2] + W_Lapl[(dimX*dimY)*k + j*dimX+i1] - 2*W_Lapl[index];
gradYYc = W_Lapl[(dimX*dimY)*k + j2*dimX+i] + W_Lapl[(dimX*dimY)*k + j1*dimX+i] - 2*W_Lapl[index];
gradZZc = W_Lapl[(dimX*dimY)*k2 + j*dimX+i] + W_Lapl[(dimX*dimY)*k1 + j*dimX+i] - 2*W_Lapl[index];

Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc + gradZZc) - (Output[index] - Input[index]));
}
}
}
return *Output;

#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,k,k1,k2,index,gradXXc,gradYYc,gradZZc)
for(k=0; k<dimZ; k++) {
/* symmetric boundary conditions */
k1 = k+1; if (k1 == dimZ) k1 = k-1;
k2 = k-1; if (k2 < 0) k2 = k+1;
for(j=0; j<dimY; j++) {
/* symmetric boundary conditions */
j1 = j+1; if (j1 == dimY) j1 = j-1;
j2 = j-1; if (j2 < 0) j2 = j+1;
for(i=0; i<dimX; i++) {
/* symmetric boundary conditions */
i1 = i+1; if (i1 == dimX) i1 = i-1;
i2 = i-1; if (i2 < 0) i2 = i+1;

index = (dimX*dimY)*k + j*dimX+i;

gradXXc = W_Lapl[(dimX*dimY)*k + j*dimX+i2] + W_Lapl[(dimX*dimY)*k + j*dimX+i1] - 2*W_Lapl[index];
gradYYc = W_Lapl[(dimX*dimY)*k + j2*dimX+i] + W_Lapl[(dimX*dimY)*k + j1*dimX+i] - 2*W_Lapl[index];
gradZZc = W_Lapl[(dimX*dimY)*k2 + j*dimX+i] + W_Lapl[(dimX*dimY)*k1 + j*dimX+i] - 2*W_Lapl[index];

Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc + gradZZc) - (Output[index] - Input[index]));
}}}
return *Output;
}
Loading

0 comments on commit 51e49a4

Please sign in to comment.