% FFTEST - Tests the Flatfield Data and generates the Flatfield
% correction image with the standard deviation of each pixel
%
% [Im , <sigma>, <info>] = fftest(ImStack, <fileOut>, <printer>)
%
% Im : Flatfield correction image, normalized around 1 to correct
% a pixel's' different quantum efficiency for incoming photons.
% sigma : Standard deviation of each pixel of Im (ideally, this reflects
% the counting statistics of each pixel, normalized to the mean
% of Im).
% info : struct-variable with additionanal information of the analysis
% of the flatfield:
% info.mean : Mean value of all flatfields
% info.low : lower count-acceptance tolerance
% info.up : upper count-acceptance tolerance
% info.tol : tolerance percentage of outliers
% info.bad : number of bad pixels
%
% ImStack : Flatfield data, generated by fillstack
% fileOut : optional output-"mat"-file, with the fields
% ff_dat.ff_mask = Im;
% ff_dat.ff_err = sigma;
% ff_dat.ff_info = info;
% default is ff_corr.mat
% printer : optional argument which specifies the printer. default is
% WBBA_005_1
%
% <argument> depicts an optional argument
%==========================================================================
% Further info for function FFTEST
% --------------------------------
% Function subroutine:
% -------------------
%
% HISTSUMMASK - gets the histmask from the summed histmasks of each and
% every ff image taken.
%
% POISSON - calculates a theoretical poisson distribution and delivers the
% value at position x
%
%==========================================================================
%
% FUNCTION: fftest.m
% ==========
%
% $Date: 2006/12/11 12:39:40 $
% $Author: pauli_s $
% $Revision: 1.2 $
% $Source: /import/cvs/X/PILATUS/X04SA/App/lib/X_PILATUS_X04SA_Matlab/fftest.m,v $
% $Tag: $
%
%
% Author(s): S. Pauli (SP)
% Co - Authors Roger Herger (RH), Christian Schlepuetz (CS)
% Address: Swiss Light Source (SLS)
% Paul Scherrer Institute
% CH - 5232 Villigen PSI
% Created: 2006/02/20
%
% Change Log:
% -----------
%
% 2006/02/20 (SP)
% - created
%
% 2006/02/22 (SP)
% - changed Input to just ImStack, instead of the same variables, as
% fillstack has
%
% 2006/03/09 (SP)
% - added printing and saving dialog and set the plots on 2 figures
%
% 2006/04/05 (SP)
% - optimized input and presentation, axis of first histogram are now
% changeable, output will be saved in a file
%
% 2006/05/06 (SP)
% - Saved file is a mat-file, mask and maskerror as struct, added printer
% input in input arguments
%
% 2006/12/07 (SP)
% - Made the graphical output nicer
% - Added a theoretical poissonian into the histogram plot
% - added output-variable info
%==========================================================================
%
%% Main function - fftest
% ======
function [Im sigma info]=fftest(ImStack,fileOut,printer)
%% check input arguments
% is there 1 input arguments?
error(nargchk(1, 3, nargin));
if (nargin < 3)
printer = 'WBBA_005_1' ;
%printer = 'WSLA_X04_2' ;
%printer = 'WLGA_227_1' ;
if (nargin < 2)
fileOut='ff_corr.mat';
end;
end;
% check if ImStack has 3 dimensions
if (ndims(ImStack) ~= 3)
error(strcat('Invalid input for ''Stack'' in function ', ...
' immerge.\nUse ''help immerge'' for further', ...
' information.\n'), ...
'')
end;
% Some variables read out of the input variables
[ydim xdim zdim] = size(ImStack);
%% check output argument
% are 1 to 3 output arguments specified?
error(nargoutchk(1, 3, nargout));
%% Calculating and plotting of a histogram
% with counting statistic of each pixel of the stack
fprintf(1,'Preparing histogram...');
set(0,'DefaultFigureWindowStyle','docked');
fig1 = figure;
% ---
% preparing of the stack and making a histogram
shstack = reshape(ImStack,1,[])'; % vector with cts per pixel
peak = mean(shstack); % mean value of shstack
xmax = 3*round(peak); % the maximal value of histogram bins
bin = 0:1:xmax;
[hdat,bdat]=hist(shstack,bin);
%----
% find the position of the histogram's maximal value
m=max(hdat(2:size(hdat,2)-1));
peak = find(hdat==m);
sigp = sqrt(peak);
axes('units','normalized','position',[0.06 0.72 0.38 0.22]);
bar(bdat,hdat);
grid on; hold on;
xlim([-2 xmax]);
xlabel('#counts');
ylabel('#pixels');
plot(1:1:xmax , ...
length(shstack)*poisson(1:1:xmax,find(hdat==max(hdat))),'r');
title({'Average counting-rate of each pixel';...
[' over ',num2str(zdim),' images']});
ylim([0 1.1*max( ...
length(shstack)*poisson(1:1:xmax,find(hdat==max(hdat))) )]);
clear shstack bdat hdat;
%% Input dialog for changing the axis
fprintf(1,'done\n');
axReply = input(['Do you want to change the ',...
'limits of the x- and y-axes? y/n [n]: '],'s');
if isempty(axReply)
axReply = 'n';
end;
ylimits = get(gca,'ylim');
if (axReply=='y')
reply = 'n';
axmin = -2;
axmax = xmax;
aymax = ylimits(2);
while (reply ~='y')
rpAxmin = ...
input(['Set lower x - limit [',int2str(axmin),']: ']);
if ~isempty(rpAxmin)
axmin = rpAxmin;
end;
rpAxmax = ...
input(['Set upper x - limit [',int2str(axmax),']: ']);
if ~isempty(rpAxmax)
axmax = rpAxmax;
end;
while (axmax < axmin)
rpAxmax = ...
input(['Set UPPER!!! x - limit [',int2str(axmax),']: ']);
if ~isempty(rpAxmax)
axmax = rpAxmax;
end;
end;
rpAymax = ...
input(['Set upper y - limit [',int2str(aymax),']: ']);
if ~isempty(rpAymax)
aymax = rpAymax;
end;
axis([axmin axmax 0 aymax]);
reply = input('Are the axis OK? y/n [y]: ', 's');
if isempty(reply)
reply = 'y';
end;
end;
end;
info.mean = peak;
%% Input dialog for the limits of the histogram
reply = 'n';
loLim = round((peak-6*sigp)/10)*10;
hiLim = round((peak+6*sigp)/10)*10;
while (reply ~='y')
rpLoLim = input(...
['Set lowest acceptable counting level [',int2str(loLim),']: ']);
if ~isempty(rpLoLim)
loLim = rpLoLim;
end;
line([loLim,loLim],[0,ylimits(2)],'Color','b');
rpHiLim = input(...
['Set highest acceptable counting level [',int2str(hiLim),']: ']);
if ~isempty(rpHiLim)
hiLim = rpHiLim;
end;
line([hiLim,hiLim],[0,ylimits(2)],'Color','b');
reply = input('Are these values OK? y/n [y]: ', 's');
if isempty(reply)
reply = 'y';
elseif (reply ~='y')
line([loLim,loLim],[0,ylimits(2)],'Color',[0.5 0.5 0.5]);
line([hiLim,hiLim],[0,ylimits(2)],'Color',[0.5 0.5 0.5]);
end;
end;
info.low = loLim;
info.up = hiLim;
clear peak;
%% Calculation of the Tolerance and plot the # of bad pix
fprintf(1,'Calculating Tolerance Graph...\n');
% Calculation of the Tail
[maskimgout, summedmask, nbp] = ...
histsummask(ImStack, [loLim hiLim], 1);
tol=zdim;
j=1;
pixct = 0;
n = 0;
while (tol>0)
[maskimgout, summedmask, nobadpix] = ...
histsummask(ImStack, [loLim hiLim], tol);
pixct(j) = nobadpix ;
n(j) = zdim-tol ;
fprintf(1,strcat(int2str(zdim-tol),', '));
tol=tol-1;
if (nobadpix ~=0)
if (tol < zdim-5) && (nobadpix < 2*nbp) && ...
( ( (oldnobadpix-nobadpix) /nobadpix) <0.25) && (tol~=0)
if (tol>20)
if (tol<zdim/2)
tol=1;
else
tol=tol-9;
end;
else
tol=1;
end;
end;
end;
j=j+1;
oldnobadpix=nobadpix;
if (tol==0)
[maskimgout, summedmask, nobadpix] = ...
histsummask(ImStack, [loLim hiLim], 0);
pixct(j) = nobadpix ;
n(j) = zdim-tol ;
fprintf(1,strcat(int2str(zdim-tol),'\n'));
end;
end;
axes('units','normalized','position',[0.56 0.72 0.38 0.22]);
plot((zdim-n)/zdim*100,pixct);
view(180,-90);
grid on;
title('Incidence of bad pixels vs. tolerance');
ylimits = get(gca,'ylim');
xlimits = get(gca,'xlim');
yl = (ylimits(2)-ylimits(1))/20;
text(xlimits(1),ylimits(2)-yl,{['total no. of pixels: ', ...
int2str(xdim*ydim)];'Definition of tolerance:'; ...
'e.g. 80% means 4 out of 5 times';...
'the pixel lies within the ';...
'acceptable counting rates'},...
'HorizontalAlignment','right', ...
'VerticalAlignment','top');
xlabel('Tolerance');
ylabel('Bad Pixels');
%% Input dialog for the limit of tolerance
reply = 'n';
tolLim = 98;
while (reply ~='y')
rpTolLim = ...
input(['Set Tolerance limit in % [',int2str(tolLim),'%]: ']);
if ~isempty(rpTolLim)
tolLim = rpTolLim;
end;
line([tolLim,tolLim],ylimits,'Color','r');
reply = input('Is that value OK? y/n [y]: ', 's');
if isempty(reply)
reply = 'y';
else
line([tolLim,tolLim],ylimits,'Color',[0.5 0.5 0.5]);
end;
end;
clear pixct;
info.tol = tolLim;
%% Calculating the mask with the user-given tolerance
fprintf(1,'Calculating the mask...');
[MaskImgOut, summedmask, nobadpix] = ...
histsummask(ImStack, [loLim hiLim], round(tolLim/100*zdim));
info.bad = nobadpix;
fprintf(1,'done\nCalculating the corrected flatfield and error...');
[Im,sigma] = makeffcorr(ImStack,MaskImgOut);
clear summedmask;
ff_dat.ff_mask = Im;
ff_dat.ff_err = sigma;
ff_dat.ff_info = info;
save(fileOut,'-struct','ff_dat');
%% Calculating the std.deviation of one pixel and plot the deviationfield
fprintf(1,'done\nCalculating z-direction standard-deviation...');
timeErr=std(ImStack,0,3);
sCop=mean2(timeErr);
ax1 = axes('units','normalized','position',[0.06 .39 0.6 .22]);
cTol=round(sqrt(sCop));
image(timeErr,'CDataMapping','scaled');
set(gca,'clim',[sCop-3*cTol sCop+3*cTol]);
axis([0 xdim 0 ydim]);
colormap(pilatus)
axis off;
title('Standard Deviation (SD) of each Pixel');
axes('units','normalized','position',[0.68 .39 0.04 .22]);
handle = gca;
colorbar(handle,'peer',ax1);
%% Find the worst pixel and a mean pixel and plot all their values
maxErr=max(max(timeErr.*MaskImgOut));
[indMaxErrX,indMaxErrY] = find(timeErr == maxErr);
[indMeanErrX,indMeanErrY] = ...
find((timeErr > sCop-0.0001) & (timeErr < sCop+0.0001) );
j=-3;
while isempty(indMeanErrX)
[indMeanErrX,indMeanErrY] = ...
find((timeErr.*MaskImgOut > sCop-10^(j)) & ...
(timeErr.*MaskImgOut < sCop+10^(j)) );
j=j+1;
end;
axes('units','normalized','position',[0.06 0.06 0.38 0.22]);
plot(squeeze(ImStack(indMaxErrX,indMaxErrY,:)),'Color','r');
xlim([1 zdim]);
ax1 = gca;
set(ax1,'XColor','r','YColor','r')
ax2 = axes('Position',get(ax1,'Position'),'XAxisLocation','top',...
'YAxisLocation','right','Color','none', ...
'XColor','b','YColor','b');
grid on;
hold on;
plot(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:)),...
'Color','b','Parent',ax2);
xlim([1 zdim]);
title('The worst pixel and a typical pixel');
%% Plotting a histogram over the std.deviation of every pixel
edges=0:1:3*sCop+1;
[hisTimeErr binTimeErr]=hist(reshape(timeErr,[],1),edges);
axes('units','normalized','position',[0.78 0.39 0.16 0.22]);
bar(binTimeErr,hisTimeErr);
xlim([-1 3*sCop+1]);
title({'Histogram of the';'SD error image (left)'});
set(gca,'YAxisLocation','right');
clear timeErr;
%% Plots a histogram of the worst and a typical pixel
loMeanAx=fix(min(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:))));
hiMeanAx=round(max(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:))));
binMeanAx=round((hiMeanAx-loMeanAx)/50);
hisMeanTick=loMeanAx:binMeanAx:hiMeanAx;
[hMeanDat nMeanBin]= ...
hist(squeeze(ImStack(indMeanErrX(1),indMeanErrY(1),:)),hisMeanTick);
axes('units','normalized','position',[0.56 0.06 0.16 0.22]);
bar(nMeanBin,hMeanDat);
axis([loMeanAx-1 hiMeanAx+1 0 max(hMeanDat)+1]);
title('Typical Pixel');
xlimits = get(gca,'xlim');
set(gca,'XTick',...
[xlimits(1) round((xlimits(1)+xlimits(2))/2) xlimits(2)]);
loMaxAx=fix(min(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:))));
hiMaxAx=round(max(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:))));
binMaxAx=round((hiMaxAx-loMaxAx)/50);
hisMaxTick=loMaxAx:binMaxAx:hiMaxAx;
[hMaxDat nMaxBin]=...
hist(squeeze(ImStack(indMaxErrX(1),indMaxErrY(1),:)),hisMaxTick);
axes('units','normalized','position',[0.78 0.06 0.16 0.22]);
bar(nMaxBin,hMaxDat);
axis([loMaxAx-1 hiMaxAx+1 0 max(hMaxDat)+1]);
title('Worst Pixel');
xlimits = get(gca,'xlim');
set(gca,'XTick',...
[xlimits(1) round((xlimits(2)+xlimits(1))/2) xlimits(2)]);
%% Printing and saving dialog
reply = 'n';
ansPrint='';
printerNam = printer;
PSnam = 'stats.eps';
while (reply ~='y')
ansPrint = input('\nDo you want to print the Figure? [y] :','s');
if isempty(ansPrint)
ansPrint = 'y';
end;
if (ansPrint=='y')
namPrint = input( ...
['What is the name of the printer? [',printerNam,'] : '],'s');
if isempty(namPrint)
namPrint=printerNam;
end;
namPrint2=strcat('-P',namPrint);
end;
ansPS = input( ...
'Do you want to save the figure as an eps-file? [y] :','s');
if isempty(ansPS)
ansPS = 'y';
end;
if (ansPS=='y')
namPS = input( ...
['How should the eps-file be called? [',PSnam,'] : '],'s');
if isempty(namPS)
namPS=PSnam;
end;
end;
reply = input('Are all inputs OK? y/n [y]: ', 's');
if isempty(reply)
reply = 'y';
end;
end;
if (ansPrint == 'y')
fprintf(1,'\nPrinting Statistics...');
print('-dpsc2',namPrint2,fig1) ;
fprintf(1,'done');
end;
if (ansPS == 'y')
fprintf(1,'\nSaving Postscript file...');
print('-depsc2',fig1,namPS);
fprintf(1,'done');
end;
if ishandle(fig1)
close(fig1);
end;
clear ansPrint ansPS namPS namPrint namPrint2;
clear bin binTimErr cTol edges hMaxDat ax1 ax2;
clear hiLim hiMaxAx hisMaxTick hiMeanAx hiMeanDat fig1 pos1;
clear hisMeanTick hisTimeErr indMaxErrX indMaxErrY;
clear indMeanErrX indMeanErrY j loLim loMaxAx loMeanAx maxErr maxPlot;
clear n nMaxBin nMeanBin sigp;
%% Calculating the mean flatfield
fprintf(1,'\nCalculating Meanvalue of each pixel...');
fig2 = figure;
sumStack = mean(ImStack,3);
meanFf = mean2(sumStack);
ssf = std2(sumStack)/meanFf;
fprintf(1,'done');
ax1 = axes('units','normalized','position',[0.28 .725 0.6 .21]);
imagesc(sumStack/meanFf);
axis off;
set(gca,'clim',[1-2*ssf 1+2*ssf]);
title('Normalized, mean flatfield image')
colormap(pilatus);
%% Calculating 4 stripes through the flatfield and plot them
uplim1 = round(ydim/2*rand(1)+ydim/2);
lolim1 = uplim1 - 10; dimlim1 = uplim1 - lolim1 + 1;
lolim2 = round(ydim/2*rand(1));
uplim2 = lolim2 + 10; dimlim2 = uplim2 - lolim2 + 1;
line([0 xdim],[uplim1 uplim1],'Color','r');
line([0 xdim],[uplim2 uplim2],'Color','b');
line([0 xdim],[lolim1 lolim1],'Color','r');
line([0 xdim],[lolim2 lolim2],'Color','b');
lelim1 = round(xdim/2*rand(1));
rilim1 = lelim1 + 10; dimlim3 = rilim1 - lelim1 + 1;
rilim2 = round(xdim/2*rand(1)+xdim/2);
lelim2 = rilim2 - 10; dimlim4 = rilim2 - lelim2 + 1;
line([lelim1 lelim1],[0 ydim],'Color','c');
line([lelim2 lelim2],[0 ydim],'Color','m');
line([rilim1 rilim1],[0 ydim],'Color','c');
line([rilim2 rilim2],[0 ydim],'Color','m');
axes('units','normalized','position',[0.9 .725 0.04 .21]);
handle = gca;
colorbar(handle,'peer',ax1);
axes('units','normalized','position',[0.28 0.575 0.6 0.125]);
plot(sum(sumStack(lolim1:uplim1,:)/meanFf/(dimlim1),1),'Color','r');
hold on;
plot(sum(sumStack(lolim2:uplim2,:)/meanFf/(dimlim2),1),'Color','b');
xlim([0 xdim]);
ylim([0.75 1.25]);
axes('units','normalized','position',[0.06 0.725 0.16 0.21]);
plot(sum(sumStack(:,lelim1:rilim1)/meanFf/(dimlim3),2),'Color','c');
hold on;
plot(sum(sumStack(:,lelim2:rilim2)/meanFf/(dimlim4),2),'Color','m');
xlim([0 ydim]);
ylim([0.75 1.25]);
view(90,90);
%% Calculating a corrected flatfield and plot it
clear ax
fprintf(1,'\nCalculating Corrected...');
corrStack = ImStack(:,:,round(zdim*(rand(1)))).*Im;
meanCf = mean2(corrStack);
fprintf(1,'done.');
ax1 = axes('units','normalized','position',[0.28 0.225 0.6 .21]);
imagesc(corrStack/meanCf);
axis off;
set(gca,'clim',[1-2*ssf 1+2*ssf]);
title('A flatfield-corrected, randomly chosen flatfield');
colormap(pilatus);
line([0 xdim],[uplim1 uplim1],'Color','r');
line([0 xdim],[uplim2 uplim2],'Color','b');
line([0 xdim],[lolim1 lolim1],'Color','r');
line([0 xdim],[lolim2 lolim2],'Color','b');
line([lelim1 lelim1],[0 ydim],'Color','c');
line([lelim2 lelim2],[0 ydim],'Color','m');
line([rilim1 rilim1],[0 ydim],'Color','c');
line([rilim2 rilim2],[0 ydim],'Color','m');
axes('units','normalized','position',[0.9 .225 0.04 .21]);
handle = gca;
colorbar(handle,'peer',ax1);
axes('units','normalized','position',[0.28 0.075 0.6 0.125]);
plot(sum(corrStack(lolim1:uplim1,:)/meanCf/(dimlim1),1),'Color','r');
hold on;
plot(sum(corrStack(lolim2:uplim2,:)/meanCf/(dimlim1),1),'Color','b');
xlim([0 xdim]);
ylim([0.75 1.25]);
axes('units','normalized','position',[0.06 0.225 0.16 0.21]);
plot(sum(corrStack(:,lelim1:rilim1)/meanCf/(dimlim3),2),'Color','c');
hold on;
plot(sum(corrStack(:,lelim2:rilim2)/meanCf/(dimlim4),2),'Color','m');
xlim([0 ydim]);
ylim([0.75 1.25]);
view(90,90);
clear ydim xdim zdim
%% Printing and saving dialog for the flatfield
reply = 'n';
ansPrint='';
printerNam = printer;
PSnam = 'flatfield.eps';
while (reply ~='y')
ansPrint = input('\nDo you want to print the Figure? [y] :','s');
if isempty(ansPrint)
ansPrint = 'y';
end;
if (ansPrint=='y')
namPrint = input(...
['What is the name of the Printer? [',printerNam,'] : '],'s');
if isempty(namPrint)
namPrint=printerNam;
end;
namPrint2=strcat('-P',namPrint);
end;
ansPS = input('Do you want to save the Figure as an eps-file? [y] :','s');
if isempty(ansPS)
ansPS = 'y';
end;
if (ansPS=='y')
namPS = input( ...
['How should the eps-file be called? [',PSnam,'] : '],'s');
if isempty(namPS)
namPS=PSnam;
end;
end;
reply = input('Are all inputs OK? y/n [y]: ', 's');
if isempty(reply)
reply = 'y';
end;
end;
if (ansPrint == 'y')
fprintf(1,'\nPrinting Flatfield...');
print('-dpsc2',namPrint2,fig2) ;
fprintf(1,'done');
end;
if (ansPS == 'y')
fprintf(1,'\nSaving Postscriptfile...');
print('-depsc2',fig2,namPS);
fprintf(1,'done');
end;
clear ansPrint;
clear ansPS;
clear namPS;
clear namPrint;
clear namPrint2;
if ishandle(fig2)
close(fig2);
end;
fprintf(1,'\n');
end % of main function
%% subfunction histsummask
% -----------
% HISTSUMMASK - gets the histmask from the summed histmasks of each and
% every ff image taken.
%
%
% [MaskImgOut, <SummedMask>, <nobadpix>]
% = histsummask(ImStackIn, thresh, tolerance, <MaskImgIn>)
%
% MaskImgOut: Maskimg based on tolerance factor operating on SummedMask.
% Maskimg is 1 for all pixels above "tolerance", 0 otherwise
% SummedMask: The sum of the n individual masks for each image. Each
% image mask contains 1s where pixels that lie within the
% limits set by vector "thresh", and zeros otherwise.
% Therefore SummedMask is an array that contains integer
% values between 0 and n (the number of ff images recorded,
% and therefore also the size of the 3rd dimension of
% ImStackIn.
% nobadpix : Returns the number of bad pixels in MaskImgOut
% ImStackIn : Input 3-D stack of images to mask
% thresh : Vector containing both the lower (low) and higher (high)
% thresholds, thresh = [low high]. Note: pixel intensities
% identically equal to either low or high are accepted.
% tolerance : integer number less than or equal to n, defining the
% minimum number of times any given pixel has to lie within
% the thresholds given by thresh. e.g., for a maximum number
% of bad pixel events of 10% and 100 images in the stack, set
% tolerance to 90.
% MaskImgIn : Optional input mask to select only a part of your
% ImStackIn, e.g. to work on a specific region of interest.
%
% <argument> depicts an optional argument
function [MaskImgOut, SummedMask, nobadpix] = ...
histsummask (ImStackIn, thresh, tolerance)
[ydim xdim zdim]=size(ImStackIn); zdim;
MaskImgOut = zeros ([ydim xdim]);
% sum the zdim histograms of ImStackIn to create SummedMask
a = histmask (ImStackIn, thresh);
SummedMask = sum(a,3);
% search for occasional hot pixels using the tolerance parameter
indices = (SummedMask >= tolerance);
MaskImgOut(indices) = 1; % sets bad pixels to zero
% make Maskimg with ones for all nonzero elements, count nonzero elements
% in Maskimg
nobadpix = ydim * xdim - nnz(MaskImgOut);
end % of histsummask
%% subfunction poisson
% -----------
% POISSON - calculates a theoretical poisson distribution and delivers the
% value at position x
%
% out = poisson(x,lambda)
%
% out : the value at position x of the poissonian
%
% x : the position of which you want to habe the value of the
% poissonian
% lambda : the mean-value of the distribution
function out = poisson(x,lambda)
if lambda > 100
out = 1./sqrt(2*pi*x).*(lambda./x).^x.*exp(x-lambda);
else
out = lambda.^x ./ gamma(x+1).* exp(-lambda);
end;
end