From 4bbcf34668434c661e718a66323798c68c058939 Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Mon, 7 Nov 2022 18:38:07 +0000 Subject: [PATCH] Broken version of hybrid_fLSQSR_TV --- MATLAB/Algorithms/hybrid_fLSQR_TV.m | 22 +-- Python/tigre/algorithms/__init__.py | 1 + .../algorithms/krylov_subspace_algorithms.py | 170 +++++++++++++++++- 3 files changed, 181 insertions(+), 12 deletions(-) diff --git a/MATLAB/Algorithms/hybrid_fLSQR_TV.m b/MATLAB/Algorithms/hybrid_fLSQR_TV.m index d0d8b5f8..e0f1b3ed 100644 --- a/MATLAB/Algorithms/hybrid_fLSQR_TV.m +++ b/MATLAB/Algorithms/hybrid_fLSQR_TV.m @@ -73,7 +73,7 @@ % Initialise matrices U = single(zeros(numel(proj), niter+1)); -V = single(zeros(prod(geo.nVoxel), niter)); % Malena: Check if prod(geo.nVoxel) is correct, I want size of object +V = single(zeros(prod(geo.nVoxel), niter)); Z = single(zeros(prod(geo.nVoxel), niter)); % Flexible basis M = zeros(niter+1,niter); % Projected matrix 1 T = zeros(niter+1); % Projected matrix 2 @@ -108,8 +108,8 @@ u = u/normr; U(:,1) = u(:); -if max(max(max(x0))) == 0 - W = ones(size(x0)); +if max(x0(:)) == 0 + W = ones(size(x0),'single'); else W = build_weights (x0); end @@ -200,7 +200,7 @@ rhsZk = [rhsk; zeros(ii,1)]; y = MZk\rhsZk; - errorL2(ii)=norm(rhsk - Mk*y); +% errorL2(ii)=norm(rhsk - Mk*y); d = Z(:,1:ii)*y; x = x0 + reshape(d,size(x0)) + xA0; @@ -215,11 +215,11 @@ qualMeasOut(:,ii)=Measure_Quality(x_prev,x,QualMeasOpts); end aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); - resL2(ii)=im3Dnorm(aux,'L2'); -% if ii>1 && resL2(ii)>resL2(ii-1) -% disp(['Algorithm stoped in iteration ', num2str(ii),' due to loss of ortogonality.']) -% return; -% end + errorL2(ii)=im3Dnorm(aux,'L2'); + if ii>1 && resL2(ii)>resL2(ii-1) + disp(['Algorithm stoped in iteration ', num2str(ii),' due to loss of ortogonality.']) + return; + end if (ii==1 && verbose) expected_time=toc*niter; disp('hybrid fLSQR TV'); @@ -372,9 +372,9 @@ function out = mvpE(k_aux, x , transp_flag) if strcmp(transp_flag,'transp') - out = x(:) - k_aux(:)*(ones(size(x(:)))'*x(:)); + out = x(:) - k_aux(:)*sum(x(:)); elseif strcmp(transp_flag,'notransp') - out = x(:) - ones(size(x(:)))*(k_aux(:)'*x(:)); + out = x(:) - (k_aux(:)'*x(:)); end end diff --git a/Python/tigre/algorithms/__init__.py b/Python/tigre/algorithms/__init__.py index 40f77c86..0fe068ca 100644 --- a/Python/tigre/algorithms/__init__.py +++ b/Python/tigre/algorithms/__init__.py @@ -15,6 +15,7 @@ from .krylov_subspace_algorithms import hybrid_lsqr from .krylov_subspace_algorithms import lsmr from .krylov_subspace_algorithms import irn_tv_cgls +from .krylov_subspace_algorithms import hybrid_flsqr_tv from .krylov_subspace_algorithms import ab_gmres from .krylov_subspace_algorithms import ba_gmres from .pocs_algorithms import asd_pocs diff --git a/Python/tigre/algorithms/krylov_subspace_algorithms.py b/Python/tigre/algorithms/krylov_subspace_algorithms.py index 9e419b10..f493e74d 100644 --- a/Python/tigre/algorithms/krylov_subspace_algorithms.py +++ b/Python/tigre/algorithms/krylov_subspace_algorithms.py @@ -3,12 +3,15 @@ import time import copy import numpy as np +import scipy as sp import tigre from tigre.algorithms.iterative_recon_alg import IterativeReconAlg from tigre.algorithms.iterative_recon_alg import decorator from tigre.utilities.Atb import Atb from tigre.utilities.Ax import Ax import tigre.algorithms as algs +import scipy.sparse.linalg + if hasattr(time, "perf_counter"): default_timer = time.perf_counter @@ -689,4 +692,169 @@ def run_main_iter(self): self.res=self.__compute_res__(self.res,w[0:-1],y) return self.res -ba_gmres = decorator(BA_GMRES, name="ba_gmres") \ No newline at end of file +ba_gmres = decorator(BA_GMRES, name="ba_gmres") + + +class hybrid_fLSQR_TV(IterativeReconAlg): + __doc__ = ( + " hybrid_fLSQR_TV solves the CBCT problem using preconditioned hybrid flexuble LSQR with TV regularization\n" + " AB_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.__U__ = np.zeros((self.niter+1,np.prod(self.geo.nDetector)*len(self.angles)),dtype=np.float32) + self.__V__ = np.zeros((self.niter,(np.prod(self.geo.nVoxel))),dtype=np.float32) + self.__Z__ = np.zeros((self.niter,(np.prod(self.geo.nVoxel))),dtype=np.float32) + + self.__M__ = np.zeros((self.niter,self.niter+1),dtype=np.float32) #% Projected matrix + self.__T__ = np.zeros((self.niter,self.niter),dtype=np.float32) #% Projected matrix + self.__proj_rhs__ = np.zeros((self.niter+1,1),dtype=np.float32) #% Projected right hand side + + + + def __build_weights__(self): + Dxx=np.copy(self.res) + Dyx=np.copy(self.res) + Dzx=np.copy(self.res) + + Dxx[0:-2,:,:]=self.res[0:-2,:,:]-self.res[1:-1,:,:] + Dyx[:,0:-2,:]=self.res[:,0:-2,:]-self.res[:,1:-1,:] + Dzx[:,:,0:-2]=self.res[:,:,0:-2]-self.res[:,:,1:-1] + + return (Dxx**2+Dyx**2+Dzx**2+1e-6)**(-1/4) + + def Lx(self,W,img): + img=np.reshape(img,self.res.shape) + Dxx=np.copy(img) + Dyx=np.copy(img) + Dzx=np.copy(img) + + Dxx[0:-2,:,:]=img[0:-2,:,:]-img[1:-1,:,:] + Dyx[:,0:-2,:]=img[:,0:-2,:]-img[:,1:-1,:] + Dzx[:,:,0:-2]=img[:,:,0:-2]-img[:,:,1:-1] + + return np.stack((W*Dxx,W*Dyx,W*Dzx),axis=0) + + def Ltx(self,W,img3): + img3 =np.reshape(img3,(3,*self.res.shape)) + Wx_1 = W * img3[0,:,:,:] + Wx_2 = W * img3[1,:,:,:] + Wx_3 = W * img3[2,:,:,:] + + DxtWx_1=Wx_1 + DytWx_2=Wx_2 + DztWx_3=Wx_3 + + DxtWx_1[1:-2,:,:]=Wx_1[1:-2,:,:]-Wx_1[0:-3,:,:] + DxtWx_1[-1,:,:]=-Wx_1[-2,:,:] + + DytWx_2[:,1:-2,:]=Wx_2[:,1:-2,:]-Wx_2[:,0:-3,:] + DytWx_2[:,-1,:]=-Wx_2[:,-2,:] + + DztWx_3[:,:,1:-2]=Wx_3[:,:,1:-2]-Wx_3[:,:,0:-3] + DztWx_3[:,:,-1]=-Wx_3[:,:,-2] + + return DxtWx_1 + DytWx_2 + DztWx_3 + + def mvT(self,k_aux,x ): + return x.ravel() - k_aux.ravel()*np.sum(x.ravel()) + def mv(self,k_aux,x ): + return x.ravel() - np.dot(k_aux.ravel(),x.ravel()) + + def run_main_iter(self): + # % Initialise matrices + + # (1) Initialize + u=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + k_aux = Ax(np.ones(self.res.shape,dtype=np.float32)/np.sqrt(np.prod(self.geo.nVoxel)), self.geo, self.angles, "Siddon", gpuids=self.gpuids) + xA0 = 1/(np.sqrt(np.prod(self.geo.nVoxel))*np.linalg.norm(k_aux.ravel(),2)**2)*np.dot(k_aux.ravel(),u.ravel()) + xA0 =np.ones(self.res.shape,dtype=np.float32)*xA0 + + k_aux = 1/(np.sqrt(np.prod(self.geo.nVoxel))*np.linalg.norm(k_aux.ravel(),2)**2)*Atb(k_aux, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + + u = u - Ax(xA0, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + normr=np.linalg.norm(u.ravel(),2) + u=u/normr + self.__U__[0]=u.ravel() + + self.__proj_rhs__[0]=normr + + if np.max(self.res)==0: + W=np.ones(self.res.shape,dtype=np.float32) + else: + W=self.__build_weights__() + + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + + for i in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(i) + v=Atb(u,self.geo,self.angles,backprojection_type="matched", gpuids=self.gpuids) + for j in range(i+1): + self.__T__[i,j] = np.dot(self.__V__[j],v.ravel()) + v = v.ravel() - self.__T__[i,j]*self.__V__[j] + v=np.reshape(v,self.res.shape) + self.__T__[i,i] = np.linalg.norm(v.ravel(),2) + v=v/self.__T__[i,i] + + z=self.mvT(k_aux,v) + + L = sp.sparse.linalg.LinearOperator((np.prod(self.res.shape),np.prod(self.res.shape)*3), matvec=lambda x: self.Ltx(W,x).ravel(),rmatvec=lambda x: self.Lx(W,x).ravel()) + aux_z=sp.sparse.linalg.lsqr(L,z.ravel(),iter_lim=50) + L = sp.sparse.linalg.LinearOperator((np.prod(self.res.shape)*3,np.prod(self.res.shape)), matvec=lambda x: self.Lx(W,x).ravel(),rmatvec=lambda x: self.Ltx(W,x).ravel()) + z=sp.sparse.linalg.lsqr(L,aux_z[0].ravel(),iter_lim=50) + z=z[0].astype(np.float32) + + z=self.mv(k_aux,z) + + self.__V__[i]=v.ravel() + self.__Z__[i]=z.ravel() + + # Update U and projected matrix M + u = Ax(np.reshape(z,self.res.shape), self.geo, self.angles, "Siddon", gpuids=self.gpuids) + for j in range(i+1): + self.__M__[i,j] = np.dot(self.__U__[j],u.ravel()) + u = u.ravel() - np.dot(self.__M__[i,j],self.__U__[j]) + self.__M__[i,i+1]=np.linalg.norm(u.ravel(),2) + u=u/self.__M__[i,i+1] + u=np.reshape(u,self.proj.shape) + self.__U__[i+1]=u.ravel() + + ### Solve the regularized projected problem + # (using the DVF of the small projected matrix) + + Mk=self.__M__[0:i+1,0:i+2] + # Prepare the projected regularization term + WZ=np.zeros((i+1,(3*np.prod(self.geo.nVoxel))),dtype=np.float32) + for j in range(i+1): + # This can be more efficiently... + # DZ can be saved and updated at each iteration + WZ[j]=self.Lx(W,self.__Z__[j]).ravel() + __,ZRk =sp.linalg.qr(np.transpose(WZ),mode='economic') + ZRksq=ZRk[0:i+1,0:i+1] + rhsk=self.__proj_rhs__[0:i+2] + MZk=np.concatenate((np.transpose(Mk),self.lmbda*ZRksq)) + rhsZk=np.concatenate((rhsk,np.zeros((i+1,1),dtype=np.float32))) + y = np.linalg.lstsq(MZk, rhsZk,rcond=None) + print(y[0]) + d=np.matmul(np.transpose(self.__Z__[0:i+1]),y[0]) + + x = self.res + np.reshape(d,self.res.shape) + xA0 + self.l2l[0, i] = np.linalg.norm((self.proj - tigre.Ax(x, self.geo, self.angles, "Siddon", gpuids=self.gpuids)).ravel(),2) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + if self.re_init_at_iteration + 1 == i or not self.restart: + print("BA-GMRES exited due to divergence at iteration "+str(i)) + return x + + return x + +hybrid_flsqr_tv = decorator(hybrid_fLSQR_TV, name="hybrid_flsqr_tv")