function [FFCoeff,ZeroCountImage]=FlatFieldCoefficients(FlatFieldImage,Mask)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   FLAT FIELD COEFFICIENTS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Very simple script to take an averaged image, and get the flat-field
% coefficients from it. This is done by finding the image mean, then
% dividing this by the individual pixel values. Any pixels with
% zero value in the original flat-field image will give Inf results after
% this step, so we set these back to zero.
%
% As well as returning the flat field coefficients, we also return the
% logical array showing the dead pixels, so that we can deal with them
% later
%
% Note that when finding the mean pixel value, we want to ignore the pixels
% which are zero, since we will have to deal with these separately. Can
% also ignore masked pixels if masking is required
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin ~= 2)
    % Set mask to empty
    Mask=ones(256,256);
end

% Find the pixels with zero counts
ZeroCountImage=(FlatFieldImage==0);

CombinedMask=Mask&(~ZeroCountImage);

MeanPixelVal=mean(FlatFieldImage(CombinedMask))
ImageStdDev=std(FlatFieldImage(CombinedMask))



% Create the flat field coefficients. Do this by dividing the mean by the
% pixel values.

MeanImage=ones(256,256).*MeanPixelVal;
FFCoeff=MeanImage./FlatFieldImage;

% Then, we need to set the appropriate pixels to one using logical
% subscripting. Note that we choose one so that any alternative processing
% done to, say, the dead columns doesn't get clobbered by the FF correction
FFCoeff(ZeroCountImage)=1;


%{ 
Example code for correcting an image
Planar_edgeraw_THL401=MeanImage('MTF_mask_fullbeam1_THL401_12keV_',0,99);
Planar_fullbeam4_THL401=MeanImage('MTF_fullbeam4_THL401_12keV_',0,99);
[Planar_FFCoeff_THL401,Planar_ZeroCountImage_THL401]=FlatFieldCoefficients(Planar_fullbeam4_THL401);
Planar_edgecorrect_THL401=Planar_edgeraw_THL401.*Planar_FFCoeff_THL401;
%}