summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
authorEdoardo Pasca <edo.paskino@gmail.com>2017-10-19 14:56:20 +0100
committerEdoardo Pasca <edo.paskino@gmail.com>2017-10-19 14:56:20 +0100
commit98443072f33a1f46eb8ea5ab27741a6400970afa (patch)
tree9cadd941d17443c0546a6b7c7ca4ddc187350bd7 /src
parentd0de394cc4d2be254fc6b1c7c89571b58f7bd30d (diff)
downloadregularization-98443072f33a1f46eb8ea5ab27741a6400970afa.tar.gz
regularization-98443072f33a1f46eb8ea5ab27741a6400970afa.tar.bz2
regularization-98443072f33a1f46eb8ea5ab27741a6400970afa.tar.xz
regularization-98443072f33a1f46eb8ea5ab27741a6400970afa.zip
further progress
Diffstat (limited to 'src')
-rw-r--r--src/Python/test_reconstructor-os.py88
1 files changed, 54 insertions, 34 deletions
diff --git a/src/Python/test_reconstructor-os.py b/src/Python/test_reconstructor-os.py
index 3ee92fa..8820db7 100644
--- a/src/Python/test_reconstructor-os.py
+++ b/src/Python/test_reconstructor-os.py
@@ -176,6 +176,23 @@ if True:
# subset loop
counterInd = 1
+ geometry_type = fistaRecon.getParameter('projector_geometry')['type']
+ if geometry_type == 'parallel' or \
+ geometry_type == 'fanflat' or \
+ geometry_type == 'fanflat_vec' :
+
+ for kkk in range(SlicesZ):
+ sino_id, sinoT[kkk] = \
+ astra.creators.create_sino3d_gpu(
+ X_t[kkk:kkk+1], proj_geomSUB, vol_geom)
+ sino_updt_Sub[kkk] = sinoT.T.copy()
+
+ else:
+ sino_id, sino_updt_Sub = \
+ astra.creators.create_sino3d_gpu(X_t, proj_geomSUB, vol_geom)
+
+ astra.matlab.data3d('delete', sino_id)
+
for ss in range(fistaRecon.getParameter('subsets')):
print ("Subset {0}".format(ss))
X_old = X.copy()
@@ -186,29 +203,29 @@ if True:
CurrSubIndices = fistaRecon.getParameter('os_indices')\
[counterInd:counterInd+numProjSub-1]
proj_geomSUB['ProjectionAngles'] = angles[CurrSubIndeces]
+
+ shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
+ shape[1] = numProjSub
+ sino_updt_Sub = numpy.zeros(shape)
+
+ if geometry_type == 'parallel' or \
+ geometry_type == 'fanflat' or \
+ geometry_type == 'fanflat_vec' :
+
+ for kkk in range(SlicesZ):
+ sino_id, sinoT = astra.creators.create_sino3d_gpu (
+ X_t[kkk:kkk+1] , proj_geomSUB, vol_geom)
+ sino_updt_Sub[kkk] = sinoT.T.copy()
+
+ else:
+ # for 3D geometry (watch the GPU memory overflow in ASTRA < 1.8)
+ sino_id, sino_updt_Sub = \
+ astra.creators.create_sino3d_gpu (X_t, proj_geomSUB, vol_geom)
+
+ astra.matlab.data3d('delete', sino_id)
+
-## if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
-## fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or \
-## fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' :
-## # if the geometry is parallel use slice-by-slice
-## # projection-backprojection routine
-## #sino_updt = zeros(size(sino),'single');
-## proj_geomT = proj_geom.copy()
-## proj_geomT['DetectorRowCount'] = 1
-## vol_geomT = vol_geom.copy()
-## vol_geomT['GridSliceCount'] = 1;
-## sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
-## for kkk in range(SlicesZ):
-## sino_id, sino_updt[kkk] = \
-## astra.creators.create_sino3d_gpu(
-## X_t[kkk:kkk+1], proj_geom, vol_geom)
-## astra.matlab.data3d('delete', sino_id)
-## else:
-## # for divergent 3D geometry (watch the GPU memory overflow in
-## # ASTRA versions < 1.8)
-## #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
-## sino_id, sino_updt = astra.creators.create_sino3d_gpu(
-## X_t, proj_geom, vol_geom)
+
## RING REMOVAL
residual = fistaRecon.residual
@@ -217,20 +234,23 @@ if True:
fistaRecon.getParameter(['ring_lambda_R_L1',
'ring_alpha' , 'weights',
'Lipschitz_constant'])
- sino_updt_Sub = numpy.zeros(numpy.shape(self.pars['input_sinogram']))
if lambdaR_L1 > 0 :
print ("ring removal")
- geometry_type = fistaRecon.getParameter('projector_geometry')['type']
- if geometry_type == 'parallel' or \
- geometry_type == 'fanflat' or \
- geometry_type == 'fanflat_vec' :
- # here
- for kkk in range(SlicesZ):
- sino_id, sinoT[kkk] = \
- astra.creators.create_sino3d_gpu(
- X_t[kkk:kkk+1], proj_geomSUB, vol_geom)
- sino_updt_Sub[kkk] = sinoT.T.copy()
- astra.matlab.data3d('delete', sino_id)
+ residualSub = numpy.zeros(shape)
+## for a chosen subset
+## for kkk = 1:numProjSub
+## indC = CurrSubIndeces(kkk);
+## residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
+## sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram
+## end
+ for kkk in range(numProjSub):
+ indC = CurrSubIndices[kkk]
+ residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \
+ (sino_updt_Sub[:,kkk,:].squeeze() - \
+ sino[:,indC,:].squeeze() - alpha_ring * r_x)
+ # filling the full sinogram
+ sino_updt_FULL[:,indC,:] = sino_updt_Sub[:,kkk,:].squeeze()
+
## % if geometry is 2D use slice-by-slice projection-backprojection routine
## for kkk = 1:SlicesZ