function channel_images = lsmread( filename )

%stack = tiffread('C:\Users\Tom\Desktop\TAXOL2\1i2\Image92.lsm');
stack = tiffread( filename ); 

% number of images in stack
n_in_stack = size( stack, 2 ); 

n_channels = length( stack(1).data); 
% stack(1).data is a cell array of images

image_size = size( stack(1).data{1} );

% now extract all images for each channel in turn 
channel_images = zeros(image_size(1),image_size(2), n_channels);
for channelx = 1:n_channels
   all_images = zeros( image_size(1), image_size(2), n_in_stack );
   for stackx = 1:n_in_stack
      all_images(:,:,stackx) = double( stack(stackx).data{channelx} );
   end
   % assumes that we want the max of all images in stack at each pixel
   channel_images(:,:,channelx) = max( all_images, [], 3 );
   
end

% sanity check
% figure
% colormap gray
% for stackx = 1:6
%     subplot(2,3,stackx)
%     imagesc( all_images(:,:,stackx) ); 
% end
% clear all_images  % to save any possible confusion and to save memory
% 
% figure
% colormap gray
% for channelx = 1:n_channels
%     subplot(1,n_channels, channelx)
%     imagesc(channel_images{channelx});
% end
