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