"smartass" projection

classic Classic list List threaded Threaded
4 messages Options
ronan chereau ronan chereau
Reply | Threaded
Open this post in threaded view
|

"smartass" projection

*****
To join, leave or search the confocal microscopy listserv, go to:
http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
*****

Hi fellows,

Following the post on smart average
(http://labrigger.com/blog/?s=smart+average), I propose another way to
perform a projection of a z-stack that does not involve any thresolding or
any filtering.
As when you select the brightest pixel of the z-serie to obtain a maximum
projection, my idea is to look for the z-plane from which the brightest
pixel was taken and then average this pixel with the ones from the plane
before and after.
As a result, stochastic noise is greatly reduced and the real signal is
consolidated knowing that it's always contained in more than one plane.

I think some could be interested in this.
regards

Ronan Chéreau

Here is the matlab code to perform the projection (2 functions: 1 to do the
projection and 1 that is used by the first one to import a imagej like tiff
stack):



function result=smartass_projection(tiffstack,name_result,bin)


%%%% Ronan Chéreau %%%%% 121030 %%%%
%Do a max projection but special.
%Look for the max pixel thru the stack + take the pix before and after in
the z-dimension and
%average them to give the final image

%INPUTS: tiffstack     name of the tiff stack
%        name_result   specify a name for the smartass projection image
%        bin           number of pixels to bin around zmax (1 makes an
average of 3 pixels)

%REQUIRES: being in the directory that contains the tiff stack
%          import_tiffstack.m

if nargin==2, bin=1; end
if nargin==3, end

nb_frames=size(tiffstack,3);

max_proj=zeros(size(tiffstack,1),size(tiffstack,2));
numero=zeros(size(tiffstack,1),size(tiffstack,2));
result=zeros(size(tiffstack,1),size(tiffstack,2));


%Does a max projection
for x=1:size(tiffstack,1);
    for y=1:size(tiffstack,2);
        max_proj(x,y)=max(tiffstack(x,y,1:nb_frames));
    end
end

%Finds the image from which the max pixel is coming from and stores it in
%numero.
for x=1:size(tiffstack,1);
    for y=1:size(tiffstack,2);
        for z=1:nb_frames;
                if tiffstack(x,y,z)==max_proj(x,y),
                    numero(x,y)=z;
                end
        end
    end
end

%Does the smartass projection: for every x,y position on the stack, it
%goes to the z-plane that has the max intensity and average it with the pixel
%from the plane before and the plane after.
for x=1:size(tiffstack,1);
    for y=1:size(tiffstack,2);
        if numero(x,y)<bin+1,
           
result(x,y)=sum(tiffstack(x,y,numero(x,y))+tiffstack(x,y,numero(x,y)+1));
        elseif numero(x,y)>nb_frames-bin,
           
result(x,y)=sum(tiffstack(x,y,numero(x,y)-1)+tiffstack(x,y,numero(x,y)));
        else
        result(x,y)=sum(tiffstack(x,y,numero(x,y)-bin:numero(x,y)+bin));
        end
    end

end

%8 or 16 bit result
maxpix=max(result(:));

if maxpix>255,
    result=uint16(result);
else
    result=uint8(result);
end

imwrite(result,name_result,'tif');






function FinalImage=import_tiffstack(FileTif)
%from "http://www.matlabtips.com/how-to-load-tiff-stacks-fast-really-fast/"

%Imports a tiff stack.
%requirement: be in the directory
InfoImage=imfinfo(FileTif);
mImage=InfoImage(1).Width;
nImage=InfoImage(1).Height;
NumberImages=length(InfoImage);
FinalImage=zeros(nImage,mImage,NumberImages,'uint16');
 
TifLink = Tiff(FileTif, 'r');
for i=1:NumberImages
   TifLink.setDirectory(i);
   FinalImage(:,:,i)=TifLink.read();
end
TifLink.close();
Zdenek Svindrych Zdenek Svindrych
Reply | Threaded
Open this post in threaded view
|

Re: "smartass" projection

*****
To join, leave or search the confocal microscopy listserv, go to:
http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
*****

Hi Ronan,

this approach can be generatised to other "convolution kernels" than ...0,
0, 1/3, 1/3, 1/3, 0, 0, ... (in z-direction, you know)

Post some links to your images so we can quickly appreciate the improvement!
I didn't try it myself, I have no confocal stacks at hand right now.

Also the piece of code you have included looks more like Visual Basic than
Matlab... You are not taking the advantage of the huge parallelisation
Matlab syntax offers. Try and rewrite it without the for-loops! I have tried
it, well, the result is not as compact as I thought it should have been.
Indeed, there are no loops and this code should run at least a million times
faster (in interpretted mode, excluding input/output times, of course). And
without the boring comments it's just seven lines of nice, simple and easily
readable code for all the confocal Matlab guys out here :-).

(PS: don't take it seriously, without a proper dataset I am not even able to
tell whether it works as intended...)

Cheers,

deden

%%%%%%% Matlab code begins here %%%%%%%%

% initialise our test image stack
src = zeros(10,10,4);

% fill it with some values
src = rand(size(src));

% find the maximum value in each XY position and the Z positions of these maxima
[maxproj, maxZpos] = max(src,[],3);

% find the positions above and below these maximas
% the edge values are adjusted (and subsequently treated twice)
% this can be changed, e.g. to an average of two vaues only...
maxZposPlus = min(maxZpos+1,size(src,3));
maxZposMinus = max(maxZpos-1,1);

% create the "above" and "below" images
% this is a bit tricky, maybe there is a better way around...

% first create a list of indices of all XY positions
[allX allY] = ind2sub(size(maxproj), 1:numel(maxproj));

% now dump the pixel values using these indiceas and the appropriate Z-coords
% to check that we did the indexing right, maxproj-maxprojCheck should
return all zeros!
maxprojCheck = reshape( src(sub2ind(size(src), allX, allY, maxZpos(:)')),
size(maxproj) ); % not needed, just a check
maxprojPlus = reshape( src(sub2ind(size(src), allX, allY, maxZposPlus(:)')),
size(maxproj) );
maxprojMinus = reshape( src(sub2ind(size(src), allX, allY,
maxZposMinus(:)')), size(maxproj) );

% and now the final image
maxproj_smartass = (maxproj + maxprojPlus + maxprojMinus) ./ 3;

%%%%%%% Matlab code ends here %%%%%%%%



On Fri, 21 Dec 2012 05:34:50 -0600, Ronan Chéreau <[hidden email]>
wrote:

>*****
>To join, leave or search the confocal microscopy listserv, go to:
>http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
>*****
>
>Hi fellows,
>
>Following the post on smart average
>(http://labrigger.com/blog/?s=smart+average), I propose another way to
>perform a projection of a z-stack that does not involve any thresolding or
>any filtering.
>As when you select the brightest pixel of the z-serie to obtain a maximum
>projection, my idea is to look for the z-plane from which the brightest
>pixel was taken and then average this pixel with the ones from the plane
>before and after.
>As a result, stochastic noise is greatly reduced and the real signal is
>consolidated knowing that it's always contained in more than one plane.
>
>I think some could be interested in this.
>regards
>
>Ronan Chéreau
>
>Here is the matlab code to perform the projection (2 functions: 1 to do the
>projection and 1 that is used by the first one to import a imagej like tiff
>stack):
>
>
George McNamara George McNamara
Reply | Threaded
Open this post in threaded view
|

Re: "smartass" projection /// PiMP your cells

In reply to this post by ronan chereau
*****
To join, leave or search the confocal microscopy listserv, go to:
http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
*****

Hi Ronan,

But if the brightest value is noise (ex. high gain PMT, no light), you
are averaging a bright noise value with two much dimmer values. You
would be better off taking the median (for 3 pixels, the middle value).

///

My University tech transfer office arranged with VIB a license for the
Munck et al 2012 PiMP plugin described at
http://jcs.biologists.org/content/early/2012/02/21/jcs.098939
(now free online access for the paper).
I've been using an early version of the plugin for months - it works.
PiMPing 25 confocal frames (no averaging for each frame) gives a much
nicer result than simply averaging the 25 timepoints - an extreme
version of what Ronan is talking about.


George
p.s. for PiMPing LSM710 and SP5 inverted data I like to use for 1.4 NA
objective lens: 30 nm pixel size (1/2 the 60 nm I usually use to meet
sampling criteria), 25 frames, no averaging, modest gain (500, 550 or
600). PiMP runs as a plugin in Fiji ImageJ. simplest to acquire and save
one time series, one channel per LSM or LIF file. The PiMP plugin
reports the amount of photobleaching - usually is less than the 5% they
estimate as being optimal. Works fine anyway.


On 12/21/2012 6:34 AM, Ronan Chéreau wrote:

> *****
> To join, leave or search the confocal microscopy listserv, go to:
> http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
> *****
>
> Hi fellows,
>
> Following the post on smart average
> (http://labrigger.com/blog/?s=smart+average), I propose another way to
> perform a projection of a z-stack that does not involve any thresolding or
> any filtering.
> As when you select the brightest pixel of the z-serie to obtain a maximum
> projection, my idea is to look for the z-plane from which the brightest
> pixel was taken and then average this pixel with the ones from the plane
> before and after.
> As a result, stochastic noise is greatly reduced and the real signal is
> consolidated knowing that it's always contained in more than one plane.
>
> I think some could be interested in this.
> regards
>
> Ronan Chéreau
>
> Here is the matlab code to perform the projection (2 functions: 1 to do the
> projection and 1 that is used by the first one to import a imagej like tiff
> stack):
>
>
>
> function result=smartass_projection(tiffstack,name_result,bin)
>
>
> %%%% Ronan Chéreau %%%%% 121030 %%%%
> %Do a max projection but special.
> %Look for the max pixel thru the stack + take the pix before and after in
> the z-dimension and
> %average them to give the final image
>
> %INPUTS: tiffstack     name of the tiff stack
> %        name_result   specify a name for the smartass projection image
> %        bin           number of pixels to bin around zmax (1 makes an
> average of 3 pixels)
>
> %REQUIRES: being in the directory that contains the tiff stack
> %          import_tiffstack.m
>
> if nargin==2, bin=1; end
> if nargin==3, end
>
> nb_frames=size(tiffstack,3);
>
> max_proj=zeros(size(tiffstack,1),size(tiffstack,2));
> numero=zeros(size(tiffstack,1),size(tiffstack,2));
> result=zeros(size(tiffstack,1),size(tiffstack,2));
>
>
> %Does a max projection
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          max_proj(x,y)=max(tiffstack(x,y,1:nb_frames));
>      end
> end
>
> %Finds the image from which the max pixel is coming from and stores it in
> %numero.
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          for z=1:nb_frames;
>                  if tiffstack(x,y,z)==max_proj(x,y),
>                      numero(x,y)=z;
>                  end
>          end
>      end
> end
>
> %Does the smartass projection: for every x,y position on the stack, it
> %goes to the z-plane that has the max intensity and average it with the pixel
> %from the plane before and the plane after.
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          if numero(x,y)<bin+1,
>
> result(x,y)=sum(tiffstack(x,y,numero(x,y))+tiffstack(x,y,numero(x,y)+1));
>          elseif numero(x,y)>nb_frames-bin,
>
> result(x,y)=sum(tiffstack(x,y,numero(x,y)-1)+tiffstack(x,y,numero(x,y)));
>          else
>          result(x,y)=sum(tiffstack(x,y,numero(x,y)-bin:numero(x,y)+bin));
>          end
>      end
>
> end
>
> %8 or 16 bit result
> maxpix=max(result(:));
>
> if maxpix>255,
>      result=uint16(result);
> else
>      result=uint8(result);
> end
>
> imwrite(result,name_result,'tif');
>
>
>
>
>
>
> function FinalImage=import_tiffstack(FileTif)
> %from "http://www.matlabtips.com/how-to-load-tiff-stacks-fast-really-fast/"
>
> %Imports a tiff stack.
> %requirement: be in the directory
> InfoImage=imfinfo(FileTif);
> mImage=InfoImage(1).Width;
> nImage=InfoImage(1).Height;
> NumberImages=length(InfoImage);
> FinalImage=zeros(nImage,mImage,NumberImages,'uint16');
>
> TifLink = Tiff(FileTif, 'r');
> for i=1:NumberImages
>     TifLink.setDirectory(i);
>     FinalImage(:,:,i)=TifLink.read();
> end
> TifLink.close();
>
>    
Guy Cox-2 Guy Cox-2
Reply | Threaded
Open this post in threaded view
|

Re: "smartass" projection /// PiMP your cells

*****
To join, leave or search the confocal microscopy listserv, go to:
http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
*****

George said: " You would be better off taking the median (for 3 pixels, the middle value)"

Something I proposed some years ago ...

G.C. Cox and Colin Sheppard, 1999  Appropriate Image Processing for Confocal Microscopy.  In: P.C. Cheng, P P Hwang, J L. Wu, G Wang & H Kim (eds) Focus on Multidimensional  Microscopy.  World Scientific Publishing, Singapore, New Jersey, London & Hong Kong.  Volume 2, pp 42-54  ISBN 981-02-3992-0

My method was to use a median filter with a cruciform structuring element (compare a pixel with its immediate neighbours, not diagonal ones).  This scales linearly with the number of dimensions, rather than as the power of the number of dimensions.  Back around 1995, when I did the work, computers weren't so powerful and this was important.  

You can see an example of this on a projection of a 3D stack in my chapter in the Pawley book.  (Given that I admit the original publication is a bit obscure).  

                                                                Guy

-----Original Message-----
From: Confocal Microscopy List [mailto:[hidden email]] On Behalf Of George McNamara
Sent: Saturday, 22 December 2012 11:03 AM
To: [hidden email]
Subject: Re: "smartass" projection /// PiMP your cells

*****
To join, leave or search the confocal microscopy listserv, go to:
http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
*****

Hi Ronan,

But if the brightest value is noise (ex. high gain PMT, no light), you are averaging a bright noise value with two much dimmer values. You would be better off taking the median (for 3 pixels, the middle value).

///

My University tech transfer office arranged with VIB a license for the Munck et al 2012 PiMP plugin described at
http://jcs.biologists.org/content/early/2012/02/21/jcs.098939
(now free online access for the paper).
I've been using an early version of the plugin for months - it works.
PiMPing 25 confocal frames (no averaging for each frame) gives a much nicer result than simply averaging the 25 timepoints - an extreme version of what Ronan is talking about.


George
p.s. for PiMPing LSM710 and SP5 inverted data I like to use for 1.4 NA objective lens: 30 nm pixel size (1/2 the 60 nm I usually use to meet sampling criteria), 25 frames, no averaging, modest gain (500, 550 or 600). PiMP runs as a plugin in Fiji ImageJ. simplest to acquire and save one time series, one channel per LSM or LIF file. The PiMP plugin reports the amount of photobleaching - usually is less than the 5% they estimate as being optimal. Works fine anyway.


On 12/21/2012 6:34 AM, Ronan Chéreau wrote:

> *****
> To join, leave or search the confocal microscopy listserv, go to:
> http://lists.umn.edu/cgi-bin/wa?A0=confocalmicroscopy
> *****
>
> Hi fellows,
>
> Following the post on smart average
> (http://labrigger.com/blog/?s=smart+average), I propose another way to
> perform a projection of a z-stack that does not involve any
> thresolding or any filtering.
> As when you select the brightest pixel of the z-serie to obtain a
> maximum projection, my idea is to look for the z-plane from which the
> brightest pixel was taken and then average this pixel with the ones
> from the plane before and after.
> As a result, stochastic noise is greatly reduced and the real signal
> is consolidated knowing that it's always contained in more than one plane.
>
> I think some could be interested in this.
> regards
>
> Ronan Chéreau
>
> Here is the matlab code to perform the projection (2 functions: 1 to
> do the projection and 1 that is used by the first one to import a
> imagej like tiff
> stack):
>
>
>
> function result=smartass_projection(tiffstack,name_result,bin)
>
>
> %%%% Ronan Chéreau %%%%% 121030 %%%%
> %Do a max projection but special.
> %Look for the max pixel thru the stack + take the pix before and after
> in the z-dimension and %average them to give the final image
>
> %INPUTS: tiffstack     name of the tiff stack
> %        name_result   specify a name for the smartass projection image
> %        bin           number of pixels to bin around zmax (1 makes an
> average of 3 pixels)
>
> %REQUIRES: being in the directory that contains the tiff stack
> %          import_tiffstack.m
>
> if nargin==2, bin=1; end
> if nargin==3, end
>
> nb_frames=size(tiffstack,3);
>
> max_proj=zeros(size(tiffstack,1),size(tiffstack,2));
> numero=zeros(size(tiffstack,1),size(tiffstack,2));
> result=zeros(size(tiffstack,1),size(tiffstack,2));
>
>
> %Does a max projection
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          max_proj(x,y)=max(tiffstack(x,y,1:nb_frames));
>      end
> end
>
> %Finds the image from which the max pixel is coming from and stores it
> in %numero.
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          for z=1:nb_frames;
>                  if tiffstack(x,y,z)==max_proj(x,y),
>                      numero(x,y)=z;
>                  end
>          end
>      end
> end
>
> %Does the smartass projection: for every x,y position on the stack, it
> %goes to the z-plane that has the max intensity and average it with
> the pixel %from the plane before and the plane after.
> for x=1:size(tiffstack,1);
>      for y=1:size(tiffstack,2);
>          if numero(x,y)<bin+1,
>
> result(x,y)=sum(tiffstack(x,y,numero(x,y))+tiffstack(x,y,numero(x,y)+1));
>          elseif numero(x,y)>nb_frames-bin,
>
> result(x,y)=sum(tiffstack(x,y,numero(x,y)-1)+tiffstack(x,y,numero(x,y)));
>          else
>          result(x,y)=sum(tiffstack(x,y,numero(x,y)-bin:numero(x,y)+bin));
>          end
>      end
>
> end
>
> %8 or 16 bit result
> maxpix=max(result(:));
>
> if maxpix>255,
>      result=uint16(result);
> else
>      result=uint8(result);
> end
>
> imwrite(result,name_result,'tif');
>
>
>
>
>
>
> function FinalImage=import_tiffstack(FileTif)
> %from "http://www.matlabtips.com/how-to-load-tiff-stacks-fast-really-fast/"
>
> %Imports a tiff stack.
> %requirement: be in the directory
> InfoImage=imfinfo(FileTif);
> mImage=InfoImage(1).Width;
> nImage=InfoImage(1).Height;
> NumberImages=length(InfoImage);
> FinalImage=zeros(nImage,mImage,NumberImages,'uint16');
>
> TifLink = Tiff(FileTif, 'r');
> for i=1:NumberImages
>     TifLink.setDirectory(i);
>     FinalImage(:,:,i)=TifLink.read();
> end
> TifLink.close();
>
>