some clearing

7 files changed, 29 insertions, 296 deletions

diff --git a/Readme.md b/Readme.md
index a530c72..268bb15 100644
--- a/Readme.md
+++ b/Readme.md
@@ -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