1*c4762a1bSJed Brown% (1) Generate Object 2*c4762a1bSJed Brown% sampleName: the type and name to specify a sample object, not case-sensitive 3*c4762a1bSJed Brown% Ny*Nx: Sample object size 4*c4762a1bSJed Brown% WGT: Sample object ground truth 5*c4762a1bSJed BrownsampleName = 'Phantom'; % choose from {'Phantom', 'Brain', ''Golosio', 'circle', 'checkboard', 'fakeRod'}; not case-sensitive 6*c4762a1bSJed BrownNx = 50; Ny = Nx; WSz = [Ny, Nx]; % XH: Ny = 100 -> 10 or 50 for debug. currently assume object shape is square not general rectangle 7*c4762a1bSJed BrownWGT = createObject(sampleName, WSz); 8*c4762a1bSJed Brown 9*c4762a1bSJed Brown% (2) Scan Object: Create Forward Model and Sinograms 10*c4762a1bSJed Brown% NTheta*NTau: Sinogram Size 11*c4762a1bSJed Brown% L: Forward Model, sparse matrix of NTheta*NTau by Ny*Nx 12*c4762a1bSJed Brown% SAll: Sinogram for all scans of NTheta*NTau*NScan, SAll(:,:,n) for the nth scan 13*c4762a1bSJed BrownNTheta = 20; % sample angle #. Use odd NOT even, for display purpose of sinagram of Phantom. As even NTheta will include theta of 90 degree where sinagram will be very bright as Phantom sample has verical bright line on left and right boundary. 14*c4762a1bSJed BrownNTau = ceil(sqrt(sum(WSz.^2))); NTau = NTau + rem(NTau-Ny,2); % number of discrete beam, enough to cover object diagonal, also use + rem(NTau-Ny,2) to make NTau the same odd/even as Ny just for perfection, so that for theta=0, we have sum(WGT, 2)' and S(1, (1:Ny)+(NTau-Ny)/2) are the same with a scale factor 15*c4762a1bSJed BrownSSz = [NTheta, NTau]; 16*c4762a1bSJed BrownL = XTM_Tensor_XH(WSz, NTheta, NTau); 17*c4762a1bSJed BrownS = reshape(L*WGT(:), NTheta, NTau); 18*c4762a1bSJed Brown 19*c4762a1bSJed Brown%% Save data in petsc binary format, b = A*x 20*c4762a1bSJed Brown% save to one file 21*c4762a1bSJed BrownPetscBinaryWrite('tomographyData_A_b_xGT', L, S(:), WGT(:), 'precision', 'float64'); 22*c4762a1bSJed Brown[A2, b2, xGT2] = PetscBinaryRead('tomographyData_A_b_xGT'); 23*c4762a1bSJed Browndifference(full(A2), full(L)); 24*c4762a1bSJed Browndifference(b2, S(:)); 25*c4762a1bSJed Browndifference(xGT2, WGT(:)); 26*c4762a1bSJed Brown% Save to separate files 27*c4762a1bSJed Brown% PetscBinaryWrite('tomographySparseMatrixA', L, 'precision', 'float64'); 28*c4762a1bSJed Brown% PetscBinaryWrite('tomographyVecXGT', WGT(:), 'precision', 'float64'); 29*c4762a1bSJed Brown% PetscBinaryWrite('tomographyVecB', S(:), 'precision', 'float64'); 30*c4762a1bSJed Brown 31*c4762a1bSJed Brown%% Compare Groundtruth vs BRGN reconstruction vs (optional TwIST result) 32*c4762a1bSJed Brown% Ground truth and model 33*c4762a1bSJed Brown[A, b, xGT] = PetscBinaryRead('tomographyData_A_b_xGT'); 34*c4762a1bSJed BrownNx = sqrt(numel(xGT)); Ny = Nx; WSz = [Ny, Nx]; 35*c4762a1bSJed BrownWGT = reshape(xGT, WSz); 36*c4762a1bSJed Brown% petsc reconstruction 37*c4762a1bSJed BrownxRecBRGN = PetscBinaryRead('tomographyResult_x'); 38*c4762a1bSJed BrownWRecBRGN = reshape(xRecBRGN, WSz); 39*c4762a1bSJed Brown% Prepare for figure 40*c4762a1bSJed BrownWAll = {WGT, WRecBRGN}; 41*c4762a1bSJed BrowntitleAll = {'Ground Truth', sprintf('Reconstruction-Tao-brgn-nonnegative,psnr=%.2fdB', psnr(WRecBRGN, WGT))}; 42*c4762a1bSJed Brown 43*c4762a1bSJed Brown% May add the Matlab reconstruction using TwIST to comparison 44*c4762a1bSJed BrownisDemoMatLabReconstruction = 1; % 1/0 45*c4762a1bSJed Brownif isDemoMatLabReconstruction 46*c4762a1bSJed Brown % Reconstruction with solver from XH, with L1/TV regularizer. 47*c4762a1bSJed Brown % Need 100/500/1000+ iteration to get good/very good/excellent result with small regularizer. 48*c4762a1bSJed Brown % choose small maxSVDofA to make sure initial step size is not too small. 1.8e-6 and 1e-6 could make big difference for n=2 case 49*c4762a1bSJed Brown regType = 'L1'; % 'TV' or 'L1' % TV is better and cleaner for phantom example 50*c4762a1bSJed Brown regWt = 1e-8*max(WGT(:)); % 1e-6 to 1e-8 both good for phantom, %1e-8*max(WGT(:)) use 1e-8 for brain, especically WGT is scaled to maximum of 1 not 40 51*c4762a1bSJed Brown maxIterA = 500; % 100 is not enough? 500 is 52*c4762a1bSJed Brown maxSVDofA = 1e-6; %svds(A, 1)*1e-4; % multiply by 1e-4 to make sure it is small enough so that first step in TwIST is long enough 53*c4762a1bSJed Brown paraTwist = {'xSz', WSz, 'regFun', regType, 'regWt', regWt, 'isNonNegative', 1, 'maxIterA', maxIterA, 'xGT', xGT, 'maxSVDofA', maxSVDofA, 'tolA', 1e-8}; 54*c4762a1bSJed Brown xRecTwist = solveTwist(b, A, paraTwist{:}); 55*c4762a1bSJed Brown WRecTwist = reshape(xRecTwist, WSz); 56*c4762a1bSJed Brown WAll = [WAll, {WRecTwist}]; 57*c4762a1bSJed Brown titleAll = [titleAll, {sprintf('Reconstruction-Matlab-Twist, psnr=%.2fdB', psnr(WRecTwist, WGT))}]; 58*c4762a1bSJed Brownend 59*c4762a1bSJed Brown% 60*c4762a1bSJed Brown% show results 61*c4762a1bSJed Brownfigure(99); clf; multAxes(@imagesc, WAll); multAxes(@axis, 'image'); linkAxesXYZLimColorView; multAxes(@colorbar); 62*c4762a1bSJed BrownmultAxes(@title, titleAll); 63*c4762a1bSJed Brown 64*c4762a1bSJed Brown 65*c4762a1bSJed Brown%% test PetscBinaryWrite() and PetscBinaryRead() 66*c4762a1bSJed BrowntestPetscBinaryWriteAndRead = 0; 67*c4762a1bSJed Brownif testPetscBinaryWriteAndRead 68*c4762a1bSJed Brown A = [0.81 0.91 0.28 0.96 0.96 69*c4762a1bSJed Brown 0.91 0.63 0.55 0.16 0.49 70*c4762a1bSJed Brown 0.13 0.10 0.96 0.97 0.80]; 71*c4762a1bSJed Brown xGT = [0;0;1;0;0]; 72*c4762a1bSJed Brown b = [0.28; 0.55; 0.96]; 73*c4762a1bSJed Brown D = [-1 1 0 0 0; 74*c4762a1bSJed Brown 0 -1 1 0 0; 75*c4762a1bSJed Brown 0 0 -1 1 0; 76*c4762a1bSJed Brown 0 0 0 -1 1]; 77*c4762a1bSJed Brown PetscBinaryWrite('cs1SparseMatrixA', A, 'precision', 'float64'); % do NOT need to convert A to sparse, always write as sparse matrix 78*c4762a1bSJed Brown [A2, b2, xGT2] = PetscBinaryRead('cs1Data_A_b_xGT'); 79*c4762a1bSJed Brownend