1% (1) Generate Object 2% sampleName: the type and name to specify a sample object, not case-sensitive 3% Ny*Nx: Sample object size 4% WGT: Sample object ground truth 5sampleName = 'Phantom'; % choose from {'Phantom', 'Brain', ''Golosio', 'circle', 'checkboard', 'fakeRod'}; not case-sensitive 6Nx = 50; Ny = Nx; WSz = [Ny, Nx]; % XH: Ny = 100 -> 10 or 50 for debug. currently assume object shape is square not general rectangle 7WGT = createObject(sampleName, WSz); 8 9% (2) Scan Object: Create Forward Model and Sinograms 10% NTheta*NTau: Sinogram Size 11% L: Forward Model, sparse matrix of NTheta*NTau by Ny*Nx 12% SAll: Sinogram for all scans of NTheta*NTau*NScan, SAll(:,:,n) for the nth scan 13NTheta = 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. 14NTau = 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 15SSz = [NTheta, NTau]; 16L = XTM_Tensor_XH(WSz, NTheta, NTau); 17S = reshape(L*WGT(:), NTheta, NTau); 18 19%% Save data in petsc binary format, b = A*x 20% save to one file 21PetscBinaryWrite('tomographyData_A_b_xGT', L, S(:), WGT(:), 'precision', 'float64'); 22[A2, b2, xGT2] = PetscBinaryRead('tomographyData_A_b_xGT'); 23difference(full(A2), full(L)); 24difference(b2, S(:)); 25difference(xGT2, WGT(:)); 26% Save to separate files 27% PetscBinaryWrite('tomographySparseMatrixA', L, 'precision', 'float64'); 28% PetscBinaryWrite('tomographyVecXGT', WGT(:), 'precision', 'float64'); 29% PetscBinaryWrite('tomographyVecB', S(:), 'precision', 'float64'); 30 31%% Compare Groundtruth vs BRGN reconstruction vs (optional TwIST result) 32% Ground truth and model 33[A, b, xGT] = PetscBinaryRead('tomographyData_A_b_xGT'); 34Nx = sqrt(numel(xGT)); Ny = Nx; WSz = [Ny, Nx]; 35WGT = reshape(xGT, WSz); 36% petsc reconstruction 37xRecBRGN = PetscBinaryRead('tomographyResult_x'); 38WRecBRGN = reshape(xRecBRGN, WSz); 39% Prepare for figure 40WAll = {WGT, WRecBRGN}; 41titleAll = {'Ground Truth', sprintf('Reconstruction-Tao-brgn-nonnegative,psnr=%.2fdB', psnr(WRecBRGN, WGT))}; 42 43% May add the Matlab reconstruction using TwIST to comparison 44isDemoMatLabReconstruction = 1; % 1/0 45if isDemoMatLabReconstruction 46 % Reconstruction with solver from XH, with L1/TV regularizer. 47 % Need 100/500/1000+ iteration to get good/very good/excellent result with small regularizer. 48 % 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 regType = 'L1'; % 'TV' or 'L1' % TV is better and cleaner for phantom example 50 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 maxIterA = 500; % 100 is not enough? 500 is 52 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 paraTwist = {'xSz', WSz, 'regFun', regType, 'regWt', regWt, 'isNonNegative', 1, 'maxIterA', maxIterA, 'xGT', xGT, 'maxSVDofA', maxSVDofA, 'tolA', 1e-8}; 54 xRecTwist = solveTwist(b, A, paraTwist{:}); 55 WRecTwist = reshape(xRecTwist, WSz); 56 WAll = [WAll, {WRecTwist}]; 57 titleAll = [titleAll, {sprintf('Reconstruction-Matlab-Twist, psnr=%.2fdB', psnr(WRecTwist, WGT))}]; 58end 59% 60% show results 61figure(99); clf; multAxes(@imagesc, WAll); multAxes(@axis, 'image'); linkAxesXYZLimColorView; multAxes(@colorbar); 62multAxes(@title, titleAll); 63 64 65%% test PetscBinaryWrite() and PetscBinaryRead() 66testPetscBinaryWriteAndRead = 0; 67if testPetscBinaryWriteAndRead 68 A = [0.81 0.91 0.28 0.96 0.96 69 0.91 0.63 0.55 0.16 0.49 70 0.13 0.10 0.96 0.97 0.80]; 71 xGT = [0;0;1;0;0]; 72 b = [0.28; 0.55; 0.96]; 73 D = [-1 1 0 0 0; 74 0 -1 1 0 0; 75 0 0 -1 1 0; 76 0 0 0 -1 1]; 77 PetscBinaryWrite('cs1SparseMatrixA', A, 'precision', 'float64'); % do NOT need to convert A to sparse, always write as sparse matrix 78 [A2, b2, xGT2] = PetscBinaryRead('cs1Data_A_b_xGT'); 79end