diff options
| -rw-r--r-- | src/Python/test_reconstructor-os.py | 66 |
1 files changed, 44 insertions, 22 deletions
diff --git a/src/Python/test_reconstructor-os.py b/src/Python/test_reconstructor-os.py index 6f3721f..3ee92fa 100644 --- a/src/Python/test_reconstructor-os.py +++ b/src/Python/test_reconstructor-os.py @@ -139,7 +139,10 @@ if True: ## additional for proj_geomSUB = proj_geom.copy() - fistaRecon.residual2 = numpy.zeros(numpy.shape(self.pars['input_sinogram'])) + fistaRecon.residual2 = numpy.zeros(numpy.shape(fistaRecon.pars['input_sinogram'])) + residual2 = fistaRecon.residual2 + sino_updt_FULL = residual.copy() + print ("starting iterations") ## % Outer FISTA iterations loop for i in range(fistaRecon.getParameter('number_of_iterations')): @@ -153,11 +156,23 @@ if True: # hence additional work is required one solution is to work with a full # sinogram at times - ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 - if (lambdaR_L1 > 0) : - sino_id2, sino_updt2 = astra.creators.create_sino3d_gpu( - X, proj_geom, vol_geom) - astra.matlab.data3d('delete', sino_id2) + + SlicesZ, anglesNumb, Detectors = \ + numpy.shape(fistaRecon.getParameter('input_sinogram')) ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 + if (i > 1 and lambdaR_L1 > 0) : + for kkk in range(anglesNumb): + + residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ + ((sino_updt_FULL[:,kkk,:]).squeeze() - \ + (sino[:,kkk,:]).squeeze() -\ + (alpha_ring * r_x) + ) + r_old = fistaRecon.r.copy() + vec = residual.sum(axis = 1) + #if SlicesZ > 1: + # vec = vec[:,1,:] # 1 or 0? + r_x = fistaRecon.r_x + fistaRecon.r = (r_x - (1./L_const) * vec).copy() # subset loop counterInd = 1 @@ -165,8 +180,7 @@ if True: print ("Subset {0}".format(ss)) X_old = X.copy() t_old = t - r_old = fistaRecon.r.copy() - + # the number of projections per subset numProjSub = fistaRecon.getParameter('os_bins')[ss] CurrSubIndices = fistaRecon.getParameter('os_indices')\ @@ -198,27 +212,35 @@ if True: ## RING REMOVAL residual = fistaRecon.residual - residual2 = fistaRecon.residual2 lambdaR_L1 , alpha_ring , weights , L_const= \ fistaRecon.getParameter(['ring_lambda_R_L1', 'ring_alpha' , 'weights', 'Lipschitz_constant']) - r_x = fistaRecon.r_x - SlicesZ, anglesNumb, Detectors = \ - numpy.shape(fistaRecon.getParameter('input_sinogram')) + sino_updt_Sub = numpy.zeros(numpy.shape(self.pars['input_sinogram'])) if lambdaR_L1 > 0 : print ("ring removal") -## % the ring removal part (Group-Huber fidelity) -## % first 2 iterations do additional work reconstructing whole dataset to ensure -## % the stablility -## if (i < 3) -## [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); -## astra_mex_data3d('delete', sino_id2); -## else -## [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); -## end + 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) + +## % if geometry is 2D use slice-by-slice projection-backprojection routine +## for kkk = 1:SlicesZ +## [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); +## sino_updt_Sub(:,:,kkk) = sinoT'; +## astra_mex_data2d('delete', sino_id); +## end + +## ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 if i < 3: pass @@ -254,7 +276,7 @@ if True: if i < 3: residual[:,kkk,:] = residual2[:,indC,:].squeeze() else: - residual(:,kkk,:) = \ + residual[:,kkk,:] = \ weights[:,indC,:].squeeze() * sino_updt[:,kkk,:].squeeze() - \ sino[:,indC,:].squeeze() - alpha_ring * fistaRecon.r_x #squeeze(weights(:,indC,:)).* \ |