diff options
| -rw-r--r-- | Readme.md | 13 | ||||
| -rw-r--r-- | demos/Demo1.m | 34 | ||||
| -rw-r--r-- | main_func/FISTA_REC.m | 33 | ||||
| -rw-r--r-- | supp/add_wedges.m | 35 | ||||
| -rw-r--r-- | supp/filtersinc.m | 28 | ||||
| -rw-r--r-- | supp/ssim_index.m | 181 | ||||
| -rw-r--r-- | supp/subplot_tight.m | 1 |
7 files changed, 29 insertions, 296 deletions
@@ -1,4 +1,4 @@ -# FISTA Reconstruction (Daniil Kazanteev) +# FISTA Reconstruction (Daniil Kazantsev) # General Description @@ -14,7 +14,6 @@ Software for reconstructing 2D/3D x-ray and neutron tomography datasets. The dat ### Demos: * Demo1: Synthetic phantom reconstruction with noise, stripes and zingers - * Demo2: Synthetic phantom reconstruction with noise, stripes, zingers, and the missing wedges * DemoRD1: Real data reconstruction from sino_basalt.mat (see Data) * DemoRD2: Real data reconstruction from sino3D_dendrites.mat (see Data) @@ -26,7 +25,7 @@ Software for reconstructing 2D/3D x-ray and neutron tomography datasets. The dat ### Main modules: * FISTA_REC.m – Matlab function to perform FISTA-based reconstruction - * FISTA_TV.c – C-omp function to solve for the weighted TV term using FISTA + * FGP_TV.c – C-omp function to solve for the weighted TV term using FGP * SplitBregman_TV.c – C-omp function to solve for the weighted TV term using Split-Bregman * LLT_model.c – C-omp function to solve for the weighted LLT [3] term using explicit scheme * studentst.m – Matlab function to calculate Students t penalty with 'auto-tuning' @@ -34,17 +33,13 @@ Software for reconstructing 2D/3D x-ray and neutron tomography datasets. The dat ### Supplementary: * zing_rings_add.m Matlab script to generate proj. data, add noise, zingers and stripes - * add_wedges.m script to add the missing wedge to existing sinogram - * my_red_yellowMAP.mat – nice colormap for the phantom + * my_red_yellowMAP.mat – nice colormap for the phantom * RMSE.m – Matlab function to calculate Root Mean Square Error - * subplot_tight – visualizing better subplots - * ssim_index – ssim calculation - + ### Practical advices: * Full 3D reconstruction provides much better results than 2D. In the case of ring artifacts, 3D is almost necessary * Depending on data it is better to use TV-LLT combination in order to achieve piecewise-smooth solution. The DemoRD2 shows one possible example when smoother surfaces required. * L (Lipshitz constant) if tweaked can lead to faster convergence than automatic values - * Convergence is normally much faster when using Fourier filtering before backprojection * Students’t penalty is generally quite stable in practice, however some tweaking of L might require for the real data * You can choose between SplitBregman-TV and FISTA-TV modules. The former is slower but requires less memory (for 3D volume U it can take up to 6 x U), the latter is faster but can take more memory (for 3D volume U it can take up to 11 x U). Also the SplitBregman is quite good in improving contrast. diff --git a/demos/Demo1.m b/demos/Demo1.m index 486b97c..3d57795 100644 --- a/demos/Demo1.m +++ b/demos/Demo1.m @@ -19,7 +19,7 @@ addpath('../main_func/'); addpath('../supp/');
load phantom_bone512.mat % load the phantom
-load my_red_yellowMAP.mat % load the colormap
+load my_red_yellowMAP.mat % load the colormap
% load sino1.mat; % load noisy sinogram
N = 512; % the size of the tomographic image NxN
@@ -67,12 +67,12 @@ error_FISTA = output.Resid_error; obj_FISTA = output.objective; figure(2); clf
%set(gcf, 'Position', get(0,'Screensize'));
-subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]); title('residual'); colorbar;
+subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar;
+subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]); title('residual'); colorbar;
colormap(cmapnew);
figure(3); clf
-subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA); title('RMSE plot'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA); title('Objective plot'); colorbar;
+subplot(1,2,1, [0.05 0.05]); plot(error_FISTA); title('RMSE plot'); colorbar;
+subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA); title('Objective plot'); colorbar;
colormap(cmapnew);
%%
fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...');
@@ -94,12 +94,12 @@ fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction is:', min(error_ error_FISTA_TV = output.Resid_error; obj_FISTA_TV = output.objective;
figure(4); clf
-subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]); title('residual'); colorbar;
+subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar;
+subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]); title('residual'); colorbar;
colormap(cmapnew);
figure(5); clf
-subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_TV); title('RMSE plot'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV); title('Objective plot'); colorbar;
+subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_TV); title('RMSE plot'); colorbar;
+subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV); title('Objective plot'); colorbar;
colormap(cmapnew);
%%
fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');
@@ -123,13 +123,13 @@ fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FI error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective;
figure(6); clf
-subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]); title('residual'); colorbar;
+subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;
+subplot(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]); title('residual'); colorbar;
colormap(cmapnew);
figure(7); clf
-subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_GH_TV); title('RMSE plot'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_GH_TV); title('Objective plot'); colorbar;
+subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_GH_TV); title('RMSE plot'); colorbar;
+subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_GH_TV); title('Objective plot'); colorbar;
colormap(cmapnew);
%%
fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');
@@ -155,13 +155,13 @@ error_FISTA_student_TV = output.Resid_error; obj_FISTA_student_TV = output.objec figure(8);
set(gcf, 'Position', get(0,'Screensize'));
-subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]); title('residual'); colorbar;
+subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;
+subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]); title('residual'); colorbar;
colormap(cmapnew);
figure(9);
-subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV); title('RMSE plot'); colorbar;
-subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV); title('Objective plot'); colorbar;
+subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV); title('RMSE plot'); colorbar;
+subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV); title('Objective plot'); colorbar;
colormap(cmapnew);
%%
% print all RMSE's
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index e21ba60..688dcc3 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -72,7 +72,7 @@ if (isfield(params,'L_const')) L_const = params.L_const; else % using Power method (PM) to establish L constant - niter = 5; % number of iteration for PM + niter = 6; % number of iteration for PM x = rand(N,N,SlicesZ); sqweight = sqrt(weights); [sino_id, y] = astra_create_sino3d_cuda(x, proj_geom, vol_geom); @@ -145,11 +145,6 @@ if (isfield(params,'fidelity')) else fidelity = 'LS'; end -if (isfield(params,'precondition')) - precondition = params.precondition; -else - precondition = 0; -end if (isfield(params,'show')) show = params.show; else @@ -166,6 +161,7 @@ else slice = 1; end if (isfield(params,'initialize')) + % a 'warm start' with SIRT method % Create a data object for the reconstruction rec_id = astra_mex_data3d('create', '-vol', vol_geom); @@ -191,7 +187,6 @@ else X = zeros(N,N,SlicesZ, 'single'); % storage for the solution end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Resid_error = zeros(iterFISTA,1); % error vector objective = zeros(iterFISTA,1); % obhective vector @@ -218,12 +213,7 @@ if (lambdaR_L1 > 0) add_ring(:,kkk,:) = squeeze(sino(:,kkk,:)) - alpha_ring.*r_x; end - residual = weights.*(sino_updt - add_ring); - - if (precondition == 1) - residual = filtersinc(residual'); % filtering residual (Fourier preconditioning) - residual = residual'; - end + residual = weights.*(sino_updt - add_ring); vec = sum(residual,2); if (SlicesZ > 1) @@ -295,13 +285,8 @@ else %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); [ff, gr] = studentst(res_vec,1); residual = reshape(gr, Detectors, anglesNumb, SlicesZ); - end - - if (precondition == 1) - residual = filtersinc(residual'); % filtering residual (Fourier preconditioning) - residual = residual'; - end - + end + [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom); X = X_t - (1/L_const).*x_temp; astra_mex_data3d('delete', sino_id); @@ -314,7 +299,7 @@ else else objective(i) = 0.5.*norm(residual(:))^2 + f_val; end - % X = SplitBregman_TV(single(X), lambdaTV, iterTV, tol); % TV-Split Bregman regularization on CPU (memory limited) + %X = SplitBregman_TV(single(X), lambdaTV, iterTV, tol); % TV-Split Bregman regularization on CPU (memory limited) elseif ((lambdaHO > 0) && (lambdaTV == 0)) % Higher Order regularization X = LLT_model(single(X), lambdaHO, tauHO, iterHO, tol, 0); % LLT higher order model @@ -338,10 +323,8 @@ else fprintf('%s %i %s %s %.4f %s %s %.4f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i)); else fprintf('%s %i %s %s %.4f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i)); - end - - end - + end + end end output.Resid_error = Resid_error; output.objective = objective; diff --git a/supp/add_wedges.m b/supp/add_wedges.m deleted file mode 100644 index 5bc215c..0000000 --- a/supp/add_wedges.m +++ /dev/null @@ -1,35 +0,0 @@ -% create a wedge mask to simulate the missing wedge - -[rows, columns] = size(sino_zing_rings); -grayImage = ones(rows, columns, 'uint8'); -xCoords = [0 360 0]; -yCoords = [35 7 7]; -mask = poly2mask(xCoords, yCoords, rows, columns); -grayImage(mask) = 0; - -xCoords = [727 360 727]; -yCoords = [35 7 7]; -mask = poly2mask(xCoords, yCoords, rows, columns); -grayImage(mask) = 0; - -xCoords = [0 360 0]; -yCoords = [145 173 173]; -mask = poly2mask(xCoords, yCoords, rows, columns); -grayImage(mask) = 0; - -xCoords = [727 360 727]; -yCoords = [145 173 173]; -mask = poly2mask(xCoords, yCoords, rows, columns); -grayImage(mask) = 0; - -grayImage(1:7,:) = 0; -grayImage(173:end,:) = 0; - -%figure; imshow(grayImage, [0 1]); -MW_sino_artifacts = sino_zing_rings.*double(grayImage); -% !!! -% note that we do not re-calculate Dweights for MW_sino_artifacts -% if one does: Dweights = Dweights.*double(grayImage); -% then the MW artifacts will be reduced substantially, -% through weighting. However we would like to explore -% the effect of the penalty instead.
\ No newline at end of file diff --git a/supp/filtersinc.m b/supp/filtersinc.m deleted file mode 100644 index 6c29c98..0000000 --- a/supp/filtersinc.m +++ /dev/null @@ -1,28 +0,0 @@ -function g = filtersinc(PR) - - -% filtersinc.m -% -% Written by Waqas Akram -% -% "a": This parameter varies the filter magnitude response. -% When "a" is very small (a<<1), the response approximates |w| -% As "a" is increased, the filter response starts to -% roll off at high frequencies. -a = 1; - -[Length, Count] = size(PR); -w = [-pi:(2*pi)/Length:pi-(2*pi)/Length]; - -rn1 = abs(2/a*sin(a.*w./2)); -rn2 = sin(a.*w./2); -rd = (a*w)./2; -r = rn1*(rn2/rd)^2; - -f = fftshift(r); -for i = 1:Count - IMG = fft(PR(:,i)); - fimg = IMG.*f'; - g(:,i) = ifft(fimg); -end -g = real(g);
\ No newline at end of file diff --git a/supp/ssim_index.m b/supp/ssim_index.m deleted file mode 100644 index 4fa7a79..0000000 --- a/supp/ssim_index.m +++ /dev/null @@ -1,181 +0,0 @@ -function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L) - -%======================================================================== -%SSIM Index, Version 1.0 -%Copyright(c) 2003 Zhou Wang -%All Rights Reserved. -% -%This is an implementation of the algorithm for calculating the -%Structural SIMilarity (SSIM) index between two images. Please refer -%to the following paper: -% -%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image -%quality assessment: From error visibility to structural similarity" -%IEEE Transactios on Image Processing, vol. 13, no. 4, pp.600-612, -%Apr. 2004. -% -%Kindly report any suggestions or corrections to zhouwang@ieee.org -% -%---------------------------------------------------------------------- -% -%Input : (1) img1: the first image being compared -% (2) img2: the second image being compared -% (3) K: constants in the SSIM index formula (see the above -% reference). defualt value: K = [0.01 0.03] -% (4) window: local window for statistics (see the above -% reference). default widnow is Gaussian given by -% window = fspecial('gaussian', 11, 1.5); -% (5) L: dynamic range of the images. default: L = 255 -% -%Output: (1) mssim: the mean SSIM index value between 2 images. -% If one of the images being compared is regarded as -% perfect quality, then mssim can be considered as the -% quality measure of the other image. -% If img1 = img2, then mssim = 1. -% (2) ssim_map: the SSIM index map of the test image. The map -% has a smaller size than the input images. The actual size: -% size(img1) - size(window) + 1. -% -%Default Usage: -% Given 2 test images img1 and img2, whose dynamic range is 0-255 -% -% [mssim ssim_map] = ssim_index(img1, img2); -% -%Advanced Usage: -% User defined parameters. For example -% -% K = [0.05 0.05]; -% window = ones(8); -% L = 100; -% [mssim ssim_map] = ssim_index(img1, img2, K, window, L); -% -%See the results: -% -% mssim %Gives the mssim value -% imshow(max(0, ssim_map).^4) %Shows the SSIM index map -% -%======================================================================== - - -if (nargin < 2 | nargin > 5) - ssim_index = -Inf; - ssim_map = -Inf; - return; -end - -if (size(img1) ~= size(img2)) - ssim_index = -Inf; - ssim_map = -Inf; - return; -end - -[M N] = size(img1); - -if (nargin == 2) - if ((M < 11) | (N < 11)) % ͼССû塣 - ssim_index = -Inf; - ssim_map = -Inf; - return - end - window = fspecial('gaussian', 11, 1.5); % һƫ1.511*11ĸ˹ͨ˲ - K(1) = 0.01; % default settings - K(2) = 0.03; % - L = 255; % -end - -if (nargin == 3) - if ((M < 11) | (N < 11)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - window = fspecial('gaussian', 11, 1.5); - L = 255; - if (length(K) == 2) - if (K(1) < 0 | K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end - -if (nargin == 4) - [H W] = size(window); - if ((H*W) < 4 | (H > M) | (W > N)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - L = 255; - if (length(K) == 2) - if (K(1) < 0 | K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end - -if (nargin == 5) - [H W] = size(window); - if ((H*W) < 4 | (H > M) | (W > N)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - if (length(K) == 2) - if (K(1) < 0 | K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end -%% -C1 = (K(1)*L)^2; % C1Lxyá -C2 = (K(2)*L)^2; % C2ԱȶCxyá -window = window/sum(sum(window)); %˲һ -img1 = double(img1); -img2 = double(img2); - -mu1 = filter2(window, img1, 'valid'); % ͼ˲ӼȨ -mu2 = filter2(window, img2, 'valid'); % ͼ˲ӼȨ - -mu1_sq = mu1.*mu1; % Uxƽֵ -mu2_sq = mu2.*mu2; % Uyƽֵ -mu1_mu2 = mu1.*mu2; % Ux*Uyֵ - -sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; % sigmax -sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; % sigmay -sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; % sigmaxy - -if (C1 > 0 & C2 > 0) - ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2)); -else - numerator1 = 2*mu1_mu2 + C1; - numerator2 = 2*sigma12 + C2; - denominator1 = mu1_sq + mu2_sq + C1; - denominator2 = sigma1_sq + sigma2_sq + C2; - ssim_map = ones(size(mu1)); - index = (denominator1.*denominator2 > 0); - ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index)); - index = (denominator1 ~= 0) & (denominator2 == 0); - ssim_map(index) = numerator1(index)./denominator1(index); -end - -mssim = mean2(ssim_map); - -return
\ No newline at end of file diff --git a/supp/subplot_tight.m b/supp/subplot_tight.m deleted file mode 100644 index 0b0cbd5..0000000 --- a/supp/subplot_tight.m +++ /dev/null @@ -1 +0,0 @@ -function vargout=subplot_tight(m, n, p, margins, varargin)
%% subplot_tight
% A subplot function substitude with margins user tunabble parameter.
%
%% Syntax
% h=subplot_tight(m, n, p);
% h=subplot_tight(m, n, p, margins);
% h=subplot_tight(m, n, p, margins, subplotArgs...);
%
%% Description
% Our goal is to grant the user the ability to define the margins between neighbouring
% subplots. Unfotrtunately Matlab subplot function lacks this functionality, and the
% margins between subplots can reach 40% of figure area, which is pretty lavish. While at
% the begining the function was implememnted as wrapper function for Matlab function
% subplot, it was modified due to axes del;etion resulting from what Matlab subplot
% detected as overlapping. Therefore, the current implmenetation makes no use of Matlab
% subplot function, using axes instead. This can be problematic, as axis and subplot
% parameters are quie different. Set isWrapper to "True" to return to wrapper mode, which
% fully supports subplot format.
%
%% Input arguments (defaults exist):
% margins- two elements vector [vertical,horizontal] defining the margins between
% neighbouring axes. Default value is 0.04
%
%% Output arguments
% same as subplot- none, or axes handle according to function call.
%
%% Issues & Comments
% - Note that if additional elements are used in order to be passed to subplot, margins
% parameter must be defined. For default margins value use empty element- [].
% -
%
%% Example
% close all;
% img=imread('peppers.png');
% figSubplotH=figure('Name', 'subplot');
% figSubplotTightH=figure('Name', 'subplot_tight');
% nElems=17;
% subplotRows=ceil(sqrt(nElems)-1);
% subplotRows=max(1, subplotRows);
% subplotCols=ceil(nElems/subplotRows);
% for iElem=1:nElems
% figure(figSubplotH);
% subplot(subplotRows, subplotCols, iElem);
% imshow(img);
% figure(figSubplotTightH);
% subplot_tight(subplotRows, subplotCols, iElem, [0.0001]);
% imshow(img);
% end
%
%% See also
% - subplot
%
%% Revision history
% First version: Nikolay S. 2011-03-29.
% Last update: Nikolay S. 2012-05-24.
%
% *List of Changes:*
% 2012-05-24
% Non wrapping mode (based on axes command) added, to deal with an issue of disappearing
% subplots occuring with massive axes.
%% Default params
isWrapper=false;
if (nargin<4) || isempty(margins)
margins=[0.04,0.04]; % default margins value- 4% of figure
end
if length(margins)==1
margins(2)=margins;
end
%note n and m are switched as Matlab indexing is column-wise, while subplot indexing is row-wise :(
[subplot_col,subplot_row]=ind2sub([n,m],p);
height=(1-(m+1)*margins(1))/m; % single subplot height
width=(1-(n+1)*margins(2))/n; % single subplot width
% note subplot suppors vector p inputs- so a merged subplot of higher dimentions will be created
subplot_cols=1+max(subplot_col)-min(subplot_col); % number of column elements in merged subplot
subplot_rows=1+max(subplot_row)-min(subplot_row); % number of row elements in merged subplot
merged_height=subplot_rows*( height+margins(1) )- margins(1); % merged subplot height
merged_width= subplot_cols*( width +margins(2) )- margins(2); % merged subplot width
merged_bottom=(m-max(subplot_row))*(height+margins(1)) +margins(1); % merged subplot bottom position
merged_left=min(subplot_col)*(width+margins(2))-width; % merged subplot left position
pos=[merged_left, merged_bottom, merged_width, merged_height];
if isWrapper
h=subplot(m, n, p, varargin{:}, 'Units', 'Normalized', 'Position', pos);
else
h=axes('Position', pos, varargin{:});
end
if nargout==1
vargout=h;
end
\ No newline at end of file |