Skip to content

Commit

Permalink
working 3D version
Browse files Browse the repository at this point in the history
  • Loading branch information
@dkazanc
dkazanc committed Nov 14, 2019
1 parent 8997d9a commit ad14f2b
Show file tree
Hide file tree
Showing 4 changed files with 28 additions and 61 deletions.
3 changes: 2 additions & 1 deletion Demos/Python/DemoFISTA_artifacts2D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,8 @@

RecFISTA_Huber_RingModel = RectoolsIR.FISTA(noisy_zing_stripe,
huber_data_threshold=4.0,
ring_model_horiz_size=9,
ring_model_horiz_size = 9,
ring_model_vert_size = 1,
iterationsFISTA = 20,
regularisation = 'ROF_TV',
regularisation_parameter = 0.01,
Expand Down
28 changes: 1 addition & 27 deletions Demos/Python/DemoFISTA_artifacts3D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,32 +112,6 @@
ObjSize = N_size, # a scalar to define reconstructed object dimensions
device = 'gpu')
#%%
# testing RING
Forwproj = Rectools.FORWPROJ(phantom_tm)
residual = Forwproj - projData3D_norm
#residual = np.swapaxes(residual,0,1)
#%%
from tomobar.supp.addmodules import RING_WEIGHTS
ring_model_horiz_size = 13
ring_model_vert_size=0
ring_model_slices_size=0

offset_rings = RING_WEIGHTS(residual, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size)

intens_max = 70
sliceSel = 120
plt.figure()
plt.subplot(131)
plt.imshow(offset_rings[:,sliceSel,:])
plt.title('2D Projection (erroneous)')
plt.subplot(132)
plt.imshow(offset_rings[sliceSel,:,:])
plt.title('Sinogram view')
plt.subplot(133)
plt.imshow(offset_rings[:,:,sliceSel])
plt.title('Tangentogram view')
plt.show()
#%%
print ("Reconstruction using FBP from tomobar")
recNumerical= Rectools.FBP(projData3D_norm) # FBP reconstruction

Expand Down Expand Up @@ -241,7 +215,7 @@
plt.title('3D Huber Rec, sagittal')
plt.show()
#%%
# Run FISTA reconstrucion algorithm with 3D regularisation
# Run FISTA reconstrucion algorithm with 3D regularisation and a better ring model
RecFISTA_HuberRING_TV = RectoolsIR.FISTA(projData3D_norm,
iterationsFISTA = 20,
huber_data_threshold = 1.0,
Expand Down
56 changes: 24 additions & 32 deletions src/Core/functions_CPU/RingWeights_core.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,12 +62,12 @@ float RingWeights_main(float *residual, float *weights, int horiz_window_halfsiz
}
else {
/****************3D INPUT ***************/
printf("%i %i %i \n", slices, anglesDim, detectorsDim);
/*printf("%i %i %i \n", slices, anglesDim, detectorsDim);*/
#pragma omp parallel for shared(residual, weights) private(k, j, i)
for(i = 0; i<anglesDim; i++) {
for(j = 0; j<detectorsDim; j++) {
for(k = 0; k<slices; k++) {
RingWeights3D(residual, weights, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize, full_window, anglesDim, detectorsDim, slices, i, j, k);
RingWeights3D(residual, weights, (long)(horiz_window_halfsize), (long)(vert_window_halfsize), (long)(slices_window_halfsize), full_window, anglesDim, detectorsDim, slices, i, j, k);
}}}
}
return *weights;
Expand All @@ -78,7 +78,7 @@ float RingWeights_main(float *residual, float *weights, int horiz_window_halfsiz
float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j)
{
float *Values_Vec;
long k, j1, index;
long k, j1, i1, m, index;
int counter, x, y, midval;

index = i*detectorsDim+j;
Expand All @@ -87,16 +87,14 @@ float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize,

/* intiate the estimation of the backround using strictly horizontal values */
counter = 0;
long l, i1;
for(k=-horiz_window_halfsize; k < horiz_window_halfsize; k++) {
for (k=-horiz_window_halfsize; k <= horiz_window_halfsize; k++) {
j1 = j + k;
for(l=-vert_window_halfsize; l < vert_window_halfsize; l++) {
i1 = i + l;
for (m = -vert_window_halfsize; m <=vert_window_halfsize; m++) {
i1 = i + m;
if ((j1 >= 0) && (j1 < detectorsDim) && (i1 >= 0) && (i1 < anglesDim)) {
Values_Vec[counter] = residual[i1*detectorsDim+j1];
}
Values_Vec[counter] = residual[i1*detectorsDim+j1]; }
else Values_Vec[counter] = residual[index];
counter ++;
counter++;
}}

/* perform sorting of the vector array */
Expand All @@ -107,14 +105,12 @@ float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize,
}
}
}

weights[index] = residual[index] - Values_Vec[midval];

free(Values_Vec);
return *weights;
}

float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k)
float RingWeights3D(float *residual, float *weights, long horiz_window_halfsize, long vert_window_halfsize, long slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k)
{
float *Values_Vec;
long k1, j1, i1, l, v, m, index;
Expand All @@ -123,34 +119,30 @@ float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize,
midval = (int)(0.5f*full_window) - 1;
index = (detectorsDim*anglesDim*k) + i*detectorsDim+j;
Values_Vec = (float*) calloc (full_window, sizeof(float));
/*
counter = 0;
for(l=-horiz_window_halfsize; l < horiz_window_halfsize; l++) {
j1 = j + l;
for(m=-vert_window_halfsize; m < vert_window_halfsize; m++) {
i1 = i + m;
for(v=-slices_window_halfsize; v < slices_window_halfsize; v++) {
k1 = k + v;
if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices) && (i1 >= 0) && (i1 < anglesDim)) {
Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i1*detectorsDim+j1];
}
else Values_Vec[counter] = residual[index];
counter ++;
}}}
*/

counter = 0;
for (l = -horiz_window_halfsize; l <=horiz_window_halfsize; l++) {
j1 = j + l;
for (m = -vert_window_halfsize; m <=vert_window_halfsize; m++) {
i1 = i + m;
for (v = -slices_window_halfsize; v <=slices_window_halfsize; v++) {
k1 = k + v;
if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices) && (i1 >= 0) && (i1 < anglesDim)) {
Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i1*detectorsDim+j1]; }
else Values_Vec[counter] = residual[index];
counter++;
}}}

/* perform sorting of the vector array */
/*
for (x = 0; x < counter-1; x++) {
for (y = 0; y < counter-x-1; y++) {
if (Values_Vec[y] > Values_Vec[y+1]) {
swap(&Values_Vec[y], &Values_Vec[y+1]);
}
}
}

weights[index] = residual[index] - Values_Vec[midval];
*/
j1 = j + 15;
if ((j1 > 0) && (j1 < detectorsDim)) weights[index] = residual[(detectorsDim*anglesDim*k) + i*detectorsDim+j1];

free(Values_Vec);
return *weights;
Expand Down
2 changes: 1 addition & 1 deletion src/Core/functions_CPU/RingWeights_core.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ float RingWeights_main(float *residual, float *weights, int horiz_window_halfsiz
/************2D functions ***********/
float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j);
/************3D functions ***********/
float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k);
float RingWeights3D(float *residual, float *weights, long horiz_window_halfsize, long vert_window_halfsize, long slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k);
#ifdef __cplusplus
}
#endif

0 comments on commit ad14f2b

Please sign in to comment.