Code covered by the BSD License  

Highlights from
FRETSCAL

image thumbnail

FRETSCAL

by

Eric Muller

 

For the FRET analysis of images of yeast cells labeled with YFP and CFP.

fret( varargin )
function fret( varargin )
%Version 4/3/09 This is the script that does all the real work.
%Everything else is simply passing a command to this function. It can be
%run from the gui, but could also be called from a command line in a batch
%mode. This function will find the peaks, and output the results in text
%files. It also cuts up the main input images into smaller pieces, each
%containing at their center one of the found peaks. These cut ups can then
%be inspected in the walk thru.

global fail error FileDirKnown % passed to search2_callback
global rings lastbring firstbring active
global posit FWHMconst biggestring % passed to the gaussian minimizer
dirsep = filesep;
fretversion = 'April 2, 2009';
fail = 1;
% operational flags
failaction    = 0;            % keep == failaction discard
writesmalltiff = 1;           % write out the small tifs
debug     = 0;
FileDirKnown = 0;
error     = 0;
datadir = 'out1';             % where to put results
spotThreshold = 400;          % background value used as a threshold to filter image and remove patches of high intensity that interfere with the detection of spots of interest.
% operation parameters
strain   = 'STRAIN';
acceptor = 'YFP_TAGGED';
donor    = 'CFP_TAGGED';
YFPspill     = .09;            % default
CFPspill     = .56;            %default
bringincr    = 5;
bringthick   = 3;
side         = 19;
maxmovement  = 3;             % allow brightest spot to move between images by this many pixels
maxpeaks     = 30;            % maximum spots per image
separation   = 10;           % minimum distance in pixels between spots(maybe odd at this time)
lowintensity = .50;           % dot of less than this fraction of the first dot intensity are discarded
% criteria for acceptance
YFPsignal2noise  = 1.30;       % dot intensity must this number of times greater than background
FRETsignal2noise = 1.00;
CFPsignal2noise  = 1.00;
YFPmonotonicity  = 0.05;       % within this fraction rings must be monotonic
FRETmonotonicity = 0.5;
CFPmonotonicity  = 0.5;
YFPFWHM         = 17;      % The full width at half maximum must be less than this
FRETFWHM        = 24;
CFPFWHM         = 24;
bitdepth = 16;
bpfilter = 1; %if there are bright patches in the primary image that will be filtered out. O means no, 1 means yes.
patchpixel = 1400; %the pixel intensity used by the patch filter.
context = 40;         %default is 1, no context to AOI, smallest image output
toobright = .9;    % What constitues a pixel intensity, as a fraction of max, above which camera is not linear and pixels will not be evaluated
k = 1;
FWHMminimum = .1;
frettype = 1; yfptype = 0; cfptype = 0;
limit1 = 1;
limit2 = 0;
if( isempty( varargin ) )
    error = 1;
end;
while(  k <= length( varargin ) && (error == 0) )
    if( debug == 1 )
        fprintf( 1, ' input %d %s%s\n', k, varargin{k}(1), varargin{k}(2) );
    end;
    if( varargin{ k }( 1 ) == '-' )
        if( varargin{ k }( 2 ) == 'b' ) % background variables
            if( varargin{ k }( 3 ) == 'i' )
                bringincr = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif( varargin{ k }( 3 ) == 't' )
                bringthick = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif( varargin{ k }( 3 ) == 'd' )
                bitdepth = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif( varargin{ k }( 3 ) == 'p' )
                bpfilter = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif( varargin{ k }( 3 ) == 'n' )
                patchpixel = str2double( varargin{ k+1 } );
                k = k + 1;
            else
                debug = 1;
            end;
        elseif( varargin{ k }( 2 ) == 'c' )
            context = str2double( varargin{ k+1 } );
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'd' )
            if( varargin{ k }( 3 ) == 'o' )             % dots
                maxpeaks = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif( varargin{ k }( 3 ) == 'i' )         % directory for data
                datadir = varargin{ k+1 };
                k = k + 1;
            else
                debug = 1;
            end;
        elseif( varargin{ k }( 2 ) == 'f' )
            FileDir = varargin{ k+1 };  % dir where files reside from command line
            cd(FileDir);
            k = k + 1;
            FileDirKnown = 1;
        elseif( varargin{ k }( 2 ) == 'P' )%capital
            finddot = str2double( varargin{ k+1 } );
            k = k + 1;
            if( finddot == 1 || finddot == 2 || finddot == 3 )
            else
                error = 2;
            end;
        elseif( varargin{ k }( 2 ) == 'h' )
            error = 3;
            fprintf( 1, 'fret.m version(%s)\n', fretversion );
            fprintf( 1, '-bi   => Background rings increment, in pixels, from outer edge of AOI. Default:(%d)\n', bringincr );
            fprintf( 1, '-bn   => Intensity value used by bright patch filter. Default:(%d)\n', patchpixel);
            fprintf( 1, '-bt   => Thickness of background ring, in pixels. Default:(%d)\n', bringthick );
            fprintf( 1, '-bp   => Bright patch filter on/off switch. Default:0\n');
            fprintf( 1, '-c    => Context determines the size of the output tifs in terms of pixels surrounding AOI + background. Default: (%d)\n', context );
            fprintf( 1, '-dots => Maximum number of AOIs to capture in a single image. Default:(%d)\n', maxpeaks );
            fprintf( 1, '-dir  => Directory to write outputs. Default:(%s)\n', datadir );
            fprintf( 1, '-f    => Input files directory where tif files are located.\n' );
            fprintf( 1, '-li   => Lowest intensity AOI to consider as a fraction of maximum AOI intensity. Default:(%5.2f)\n', lowintensity );
            fprintf( 1, '-ls   => Limit the search stringency limits to the primary channel. Default:(%5.2f)\n', limit1 );
            fprintf( 1, '-lc   => Limit the cull stringency limits to the primary channel Default:(%5.2f)\n', limit2 );
            fprintf( 1, '-A    => String for YFP tagged protein\n' );
            fprintf( 1, '-D    => String for CFP tagged protein\n' );
            fprintf( 1, '-S    => String for strain\n' );
            fprintf( 1, '-[C,F,Y]   => For the following flags, C, F, and Y represent CFP, FRET or YFP channel setting. Flags are case sensitive.\n' );
            fprintf( 1, '-[C,Y]s    => Spillover factors (eg. -Cs is CFP and -Ys is YFP spillover factor for FretR determination )\n' );
            fprintf( 1, '-[C,F,Y]d  => Target FWHM decay constants\n');
            fprintf( 1, '-[C,F,Y]m  => Target montonicity tolerances. \n' );
            fprintf( 1, '-[C,F,Y]n  => Target (average siginal)/(average background ratios.) \n');
            fprintf( 1, '-[C,F,Y]x  => Identifiers for image files. \n');
            fprintf( 1, '-X  => Identifiers for DIC image files. \n');
            fprintf( 1, '-E  => Experiment type: C = CFP spillover, F = FRET, or Y = YFP spillover \n');
            fprintf( 1, '-p  => Side (in pixels) of AOI\n' );
            fprintf( 1, '-P  => Set the channel to find peaks, 1=YFP, 2=FRET,3=CFP\n' );
            fprintf( 1, '-s  => Separation between AOIs. Default(%5.2f)\n', separation );
            fprintf( 1, '-T  => Areas greater than this fraction of the max intensity value of the image are not evaluated for peaks: Default: (%5.2f)\n', toobright);
            fprintf( 1, '-m  => Maximum movement of peak between images\n' );
            fprintf( 1, '-d  => Debug mode(default is off)\n' );
        elseif( varargin{ k }( 2 ) == 'l' )              %b lowest intensity for subsequent spots after the brightest is found
            if varargin{ k }( 3 ) == 's'
                limit1 = str2double( varargin{ k+1 } );
                k = k + 1;
            elseif varargin{ k }( 3 ) == 'c'
                limit2 = varargin{ k+1 };
                k = k + 1;
            elseif varargin{k}(3) == 'i'
                lowintensity = str2double( varargin{ k+1 } );
                k = k + 1;
            end
        elseif( varargin{ k }( 2 ) == 'p' )            % Where the AOI shape gets read into fret.m and name changes to side
            side = str2double( varargin{ k+1 } );
            k = k + 1;
            tside = side +1;
            if rem(tside,2)~= 0
                error = 4;
            end;
        elseif( varargin{ k }( 2 ) == 's' )            % separation between spots
            separation = str2double( varargin{ k+1 } );
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'T' )            % threshold for areas that are too bright, as a fraction of camera max
            toobright = str2double( varargin{ k+1 } );
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 't' )            % threshold for filter to remove patches of bright pixels
            spotThreshold = str2double( varargin{ k+1 } );
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'A' )            % Acceptor
            acceptor = varargin{ k+1 };
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'D' )            % Donor
            donor = varargin{ k+1 };
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'S' )            % Strain
            strain = varargin{ k+1 };
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'C' )            % CFP parameters
            if( varargin{ k }( 3 ) == 's' )
                CFPspill = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'd' )
                CFPFWHM = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'm' )
                CFPmonotonicity = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'n' )
                CFPsignal2noise = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'x' )
                Csuffix = varargin{ k+1 };
            else
                error = 5;
                fprintf( 1, 'Unknown flag command option for CFP parameters:%s\n', varargin{ k } );
            end;
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'F' )            %b FRET parameters
            if( length( varargin{ k } ) > 2 )
                if( varargin{ k }( 3 ) == 'd' )
                    FRETFWHM = str2double( varargin{ k+1 } );
                elseif( varargin{ k }( 3 ) == 'm' )
                    FRETmonotonicity = str2double( varargin{ k+1 } );
                elseif( varargin{ k }( 3 ) == 'n' )
                    FRETsignal2noise = str2double( varargin{ k+1 } );
                elseif( varargin{ k }( 3 ) == 'x' )
                    Fsuffix = varargin{ k+1 };
                else
                    fprintf( 1, 'Unknown flag command option for FRET parameters:%s\n', varargin{ k } );
                    error = 6;
                end;
            else
                fprintf( 1, 'Unknown flag command option for FRET parameters:%s\n', varargin{ k } );
                error = 7;
            end;
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'Y' )            %b YFP parameters
            if( varargin{ k }( 3 ) == 's' )
                YFPspill = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'd' )
                YFPFWHM = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'm' )
                YFPmonotonicity = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'n' )
                YFPsignal2noise = str2double( varargin{ k+1 } );
            elseif( varargin{ k }( 3 ) == 'x' )
                Ysuffix = varargin{ k+1 };
            else
                fprintf( 1, 'Unknown flag command option for YFP parameters:%s\n', varargin{ k } );
                error = 8;
            end;
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'X' )
            Dsuffix = varargin{ k+1 };
            k = k +1;
        elseif( varargin{ k }( 2 ) == 'm' )
            maxmovement = str2double(varargin{ k+1 });
            k = k + 1;
        elseif( varargin{ k }( 2 ) == 'E' )
            type = varargin{ k + 1};
            if strcmp(type,'Y')
                yfptype = 1 ;
                cfptype = 0;
                frettype = 0;
            elseif strcmp(type,'F')
                yfptype = 0 ;
                cfptype = 0;
                frettype = 1;
            elseif strcmp(type,'C')
                yfptype = 0 ;
                cfptype = 1;
                frettype = 0;
            else
                fprintf( 1, 'Unknown flag command option for Experiment parameters:%s\n', varargin{ k } );
                error = 13;
            end
            k = k + 1;
        else
            fprintf( 1, 'Unknown flag command option:%s\n', varargin{ k } );
            error = 9;
        end;
        k = k + 1;
    else
        fprintf( 1, 'Unknown command parameter:%s\n', varargin{ k } );
        k = k + 1;
        error = 10;
    end;
end;
active( 1 ) = 1; active( 2 ) = 1; active( 3 ) = 1; active(4) = 1;
%Experiment type over writes any YFPspill or CFPspill values entered with flags
if( yfptype == 1 )        % YFP spillover expt.
    active( 3 ) = 0; %CFP is not an active channel
    YFPspill = 1; %spillover factor will be calculated
    CFPspill = 0; %not used
    if( finddot == 3 )
        fprintf( 1, 'YFP spillover experiment can''t have CFP as finddot channel\n' );
        error = 11;
    end;
elseif( cfptype == 1 )    % CFP spillover expt.
    active( 1 ) = 0;  %YFP channel is not active
    CFPspill = 1; %CFP spillover factor will be calculated
    YFPspill = 0; %YFP channel not used
    if( finddot == 1 )
        fprintf( 1, 'CFP spillover experiment can''t have YFP as finddot channel\n' );
        error = 12;
    end;
else                         % FRET
end;
% Some early errors if error ~= 0 && FileDirKnown ~= 1 )
if( error == 0 )
    if( FileDirKnown == 0 )
        fprintf( 1, 'Error unknown input file\n' );
    end
elseif( error == 1 )
    fprintf( 1, 'No arguments for fret m\n' );
elseif( error == 2 )
    fprintf( 1, 'Got an illegal channel to find peaks on\n' );
elseif( error == 3 )
    %requested help summary for flags
elseif( error == 4 )
    fprintf( 1, 'The side(%d) for AOI shape must be an odd integer. \n', side );
end;
if ( error == 0 && FileDirKnown == 1)
    fdtest = isdir(FileDir);
    if(fdtest == 0);
        error = 13;
        fprintf(1, 'The directory for the files does not exist.\n');
    end;
end;
if (bringincr < 1) %there must be at least one ring of background
    bringincr =1;
end;
if (bringthick < 1)%there must be at least one ring of background
    bringthick =1;
end;
if bitdepth == 8 || bitdepth == 12 || bitdepth == 16 %bitdepth can only be 8, 12 or 16. Set to 16 by default.
else
    bitdepth = 16;
end;
maxbits = 2^(bitdepth);
toobrightmaxbits = round( toobright * maxbits );
AOIring = (side +1)/2;
rindex = AOIring-1;%Size of AOI in rings with center "ring"(pixel) is 0. Note: matlab can not make a matrix call to an index of 0, so programming takes that into consideration.
biggestring =  rindex + bringincr + bringthick -1; % Size of AOI region in rings. Center is ring 0. Background starts at the increment ring, then has thickness.
frameedge    = (biggestring) + (2 * maxmovement) + round( separation/2 ) + (side+1)/2; % Used to pad the frame of the image so we do not wander off the image matrix during analysis
lastbring    = biggestring;   % to use as background outer ring
firstbring   =  rindex + bringincr;    % to use as background inner ring
if( error == 0 && FileDirKnown == 1 ) % NO ERRORS, WE HAVE a DIRECTORY FOR INPUT, LET US SEARCH!
    point = zeros( 10, 1 ); %point in the processing of this script
    point( 1 ) = 1;
    missingfile = 0;
    if( debug == 1 )
        fprintf( 1, 'got past point( 1 ), datadir, %s;FileDir, %s \n', datadir, FileDir );
    end;
    [ s, mess, messid ] = mkdir( fullfile(FileDir,datadir ));
    if strcmp(messid,'MATLAB:MKDIR:DirectoryExists') == 1
        fprintf(1,'Removing tifs from previous analysis.\n');
        delete(strcat(FileDir,dirsep,datadir,dirsep,'*.tif')); % we are going to purge all  tif contents of output to speed up compute
    end
    if( debug ==  1 )
        fprintf( 1, 'got to first fopen %s\n', FileDir );
    end;
    fprintf(1,'Start set up.\n');
    if exist( strcat( FileDir, dirsep,datadir,dirsep, 'jmp.txt' ),'file' ) > 0
        [clean, message] = fopen( strcat( FileDir, dirsep,datadir,dirsep, 'jmp.txt' ), 'w+' ); % we are going to open jmp.txt. If a file already exists, we are going to empty the contents.
        fclose(clean);
    end
    [ fjmp, message ] = fopen( strcat( datadir, dirsep, 'jmp.txt' ), 'a+' );            % append records
    if( debug == 1 )
        fprintf( 1, 'got past second fopen %s\n', FileDir );
    end;
    fposition = ftell( fjmp );
    if( fposition == 0 )
        if frettype == 1
            fprintf( fjmp, 'Tiff_name             Strain    CFP_tag    YFP_tag    Spot YFPobs  YFPbkg   FRETobs   FRETbkg   CFPobs   CFPbkg   YFP      FRET     CFP      SPILLtot    FRETR  YFPFWHM  YFPobs/bkg  YFPMinMonTol  YFP_X  YFP_Y  FRETFWHM  FRETobs/bkg  FRETMinMonTol  FRET_X   FRET_Y  CFPFWHM CFPobs/bkg  CFPMinMonTol   CFP_X   CFP_Y\n' );
            fprintf( fjmp, 'ExptType FRET ');
            fprintf( fjmp, ' CFPspillfactor %6.3f',          CFPspill );
            fprintf( fjmp, ' YFPspillfactor %6.3f',         YFPspill );
        elseif yfptype == 1
            fprintf( fjmp, 'Tiff_name             Strain    CFP_tag    YFP_tag    Spot YFPobs  YFPbkg   FRETobs   FRETbkg   CFPobs   CFPbkg   YFP      FRET     CFP      SPILLtot    YFPspillfactor  YFPFWHM  YFPobs/bkg  YFPMinMonTol  YFP_X  YFP_Y  FRETFWHM  FRETobs/bkg  FRETMinMonTol  FRET_X   FRET_Y  CFPFWHM CFPobs/bkg  CFPMinMonTol   CFP_X   CFP_Y\n' );
            fprintf( fjmp, 'ExptType calcYFPspill ');
            fprintf( fjmp, ' CFPspillfactor %6.3f', 0.00);
            fprintf( fjmp, ' YFPspillfactor %6.3f',  1.00);
        elseif cfptype == 1
            fprintf( fjmp, 'Tiff_name             Strain    CFP_tag    YFP_tag    Spot YFPobs  YFPbkg   FRETobs   FRETbkg   CFPobs   CFPbkg   YFP      FRET     CFP      SPILLtot    CFPspillfactor  YFPFWHM  YFPobs/bkg  YFPMinMonTol  YFP_X  YFP_Y  FRETFWHM  FRETobs/bkg  FRETMinMonTol  FRET_X   FRET_Y  CFPFWHM CFPobs/bkg  CFPMinMonTol   CFP_X   CFP_Y\n' );
            fprintf( fjmp, 'ExptType calcCFPspill ');
            fprintf( fjmp, ' CFPspillfactor %6.3f', 1.00);
            fprintf( fjmp, ' YFPspillfactor %6.3f', 0.00);
        end
        fprintf( fjmp, ' PrimaryChannel %d',        finddot );
        fprintf( fjmp, ' YFPtargetFWHM %6.3f',   YFPFWHM );
        fprintf( fjmp, ' YFPtargetmono %6.3f',    YFPmonotonicity );
        fprintf( fjmp, ' YFPsignal2noise %6.3f',  YFPsignal2noise );
        fprintf( fjmp, ' FRETtargetFWHM %6.3f',  FRETFWHM );
        fprintf( fjmp, ' FRETtargetmono %6.3f',   FRETmonotonicity );
        fprintf( fjmp, ' FRETsignal2noise %6.3f', FRETsignal2noise );
        fprintf( fjmp, ' CFPtargetFWHM %6.3f',   CFPFWHM );
        fprintf( fjmp, ' CFPtargetmono %6.3f',    CFPmonotonicity );
        fprintf( fjmp, ' CFPsignal2noise %6.3f',  CFPsignal2noise );
        fprintf( fjmp, ' Bkg_incr %2d',      bringincr );
        fprintf( fjmp, ' Bkg_thick %2d',      bringthick );
        fprintf( fjmp, ' AOIshape %2d',           side );
        fprintf( fjmp, ' spotlimit %5.2f',        lowintensity );
        fprintf( fjmp, ' separation %5.2f',     separation );
        fprintf( fjmp, ' context %5.2f',     context );
        fprintf( fjmp, ' toobright %5.2f',     toobright );
        fprintf( fjmp, ' bpfilter %d',    bpfilter );
        fprintf( fjmp, ' patchpixel %d',     patchpixel );
        fprintf(fjmp,  ' maxmovement %d',  maxmovement);
        fprintf(fjmp,  ' searchlimited %d',  limit1);
        fprintf(fjmp,  ' culllimitited %d\n',  limit2);
        
    else
    end;
    filelist = dir(FileDir);
    l = length(filelist);
    Fsuffixmatch = strcat('.*',Fsuffix,'.*');
    lFsuffix = length(Fsuffix);
    filesnum = 0;
    expnum = 0;
    for j= 1:l
        if isempty(regexp(filelist(j).name,Fsuffixmatch,'match')) == 0
            filesnum = filesnum +1;
        end
    end
    Q = cell(5,filesnum);
    for j = 1:l
        if isempty(regexp(filelist(j).name,Fsuffixmatch,'match')) == 0
            expnum = expnum +1;
            Fstart = findstr(filelist(j).name,Fsuffix); %returns starting index of Fsuffix
            Fimage = filelist(j).name;
            Fleftbase = filelist(j).name(1:Fstart -1);
            Frightbase = filelist(j).name((Fstart + lFsuffix): end);
            Yimage = strcat(Fleftbase,Ysuffix,Frightbase);
            Cimage = strcat(Fleftbase,Csuffix,Frightbase);
            Dimage = strcat(Fleftbase,Dsuffix,Frightbase);
            expname = strcat(Fleftbase,'set',Frightbase);
            Q(1:5,expnum) = {Yimage;Fimage;Cimage;Dimage;expname};
        end
    end
    cd (FileDir);
    for e = 1:expnum
        Image_Set = e % just used to make following command line during processing easier
        fprintf(1,'Set up is complete. Start the analysis of %s.\nStart clock.\n',Q{5,e});
        tic
        point( 2 ) = 2;
        if active(1) == 1
            testY = exist(Q{1,e},'file');
        end
        if active(2) == 1
            testF = exist(Q{2,e},'file');
        end
        if active(3) == 1
            testC = exist(Q{3,e},'file');
        end
        if active(4) == 1
            testD = exist(Q{4,e},'file');
        end
        if active(1) == 1
            if testY~=2 || testF~=2 ||  testD ~= 2
                missingfile = 1;
                fprintf( 1, '\n ERROR! One of the image files is missing. Will skip set.\n' );
            end
        end
        if active(3) == 1
            if testC~=2 || testF~=2 ||  testD ~= 2
                fprintf( 1, '\n ERROR! One of the image files is missing. Will skip set.\n' );
                missingfile = 1;
            end
        end
        if missingfile ~= 1;
            skipall = 0;
            fprintf( 1, 'Reading in images.\n' );
            F = imread(Q{2,e});
            D = imread(Q{4,e});
            [Height0,Width0,Cell] = size( F ); %Always have FRET image which is same size as other images.
            ptest = zeros([Height0,Width0],'uint16'); %ptest is used in the filters
            if Height0 -2*frameedge <= side  || Width0 -2*frameedge <= side
                fprintf( 1, 'The analysis images by FretSCal avoids the edges of images.\n');
                fprintf(1, 'In this image set the AOI shape,%d x %d, is too large compared to the image size,%d, minus the imposed frameedge, %d. \n',side,side,Height0,frameedge);
                skipall = 1;
                fail = 1;
            end
            if( active( 3 ) == 1 )
                C = imread(Q{3,e});
            else
                C = zeros( Height0, Width0 );
            end;
            if( active( 1 ) == 1 )
                Y = imread(Q{1,e});
            else
                Y = zeros( Height0, Width0 );
            end;
            point( 3 ) = 3;
            if( debug == 1 )
                fprintf( 1, 'got past point %d\n', point(3) )
                fprintf( 1, '%s\n', FileDir );
            end;
            MaxY = max(Y(:));
            MaxF = max(F(:));
            MaxC = max(C(:));
            MaxD = max(D(:));
            MinY = min(Y(:));
            MinF = min(F(:));
            MinC = min(C(:));
            MinD = min(D(:));
            MeanY = mean(Y(:));
            MeanF = mean(F(:));
            MeanC = mean(C(:));
            fprintf(1,'Finished input of image files and got characteristics.\nStart checking simple problems.\n');
            if( ( MaxY == 1 && active(1) == 1 ) || ( MaxF == 1 && active(2) == 1 ) || ( MaxC == 1 && active(3) == 1 ) )
                fprintf( 1, 'Blank tif. Skipping %s\n',FileDir );
                skipall = 1;
                fail = 1;
            end;
            if( ( MaxY == MinY && active(1) == 1 ) || ( MaxF == MinF && active(2) == 1 ) || ( MaxC == MinC && active(3) == 1 ) )
                fprintf( 1, 'Constant tif. Skipping set\n' );
                skipall = 1;
                fail = 1;
            end;
            %             if( ( MaxY < MinY + 50 && active(1) == 1 ) || ( MaxF < MinF + 50 && active(2) == 1 ) || ( MaxC < MinC + 50 && active(3) == 1 ) )
            %                 fprintf( 1, 'Note: Max signal very low.\n' );
            %                 skipall = 0;
            %                 fail = 0;
            %             end;
            if( ( MeanY > .95*maxbits && active(1) == 1 ) || ( MeanF >.95*maxbits  && active(2) == 1 ) || ( MeanC > .95*maxbits && active(3) == 1 ) )
                fprintf( 1, 'Image is blown out. Skipping set\n' );
                skipall = 1;
                fail = 1;
            end;
            if( ( MaxY > maxbits && active(1) == 1 ) || ( MaxF >maxbits  && active(2) == 1 ) || ( MaxC > maxbits && active(3) == 1 ) || ( MaxD > maxbits && active(4) == 1 ))
                fprintf( 1, 'The max intensity of images is greater than the max for the bit depth. Check bit depth setting. Skipping set\n' );
                skipall = 1;
                fail = 1;
            end;
            if( skipall == 0 )%no errors with file set, so will first filter files, if required, and then find AOIs
                AOIregion = (biggestring*2) + 1; % the side of the AOI plus background
                dAr = AOIregion * AOIregion;
                if active(1) == 1
                    if MaxY > toobrightmaxbits
                        fprintf(1,'Masking out maxed out areas in YFP image...\n');
                        T = ordfilt2(Y, dAr, ones(AOIregion));
                        tf = T > toobrightmaxbits;
                        Y(tf) = ptest(tf);%only replace with 0 if tf is true
                        F(tf) = ptest(tf);% and we zero out the other images in the same region too.
                        if active(3) == 1
                            C(tf) = ptest(tf);
                        end
                    end
                end
                if active(2) == 1
                    if MaxF > toobrightmaxbits
                        fprintf(1,'Masking out maxed out areas in FRET image...\n');
                        T = ordfilt2(F, dAr, ones(AOIregion));
                        tf = T > toobrightmaxbits;
                        F(tf) = ptest(tf);
                        if active(3) == 1
                            C(tf) = ptest(tf);
                        end
                        if active(1) == 1
                            Y(tf) = ptest(tf);
                        end
                    end
                end
                if active(3) == 1
                    if MaxC > toobrightmaxbits
                        fprintf(1,'Masking out maxed out areas in CFP...\n');
                        T = ordfilt2(C, dAr, ones(AOIregion));
                        tf = T > toobrightmaxbits;
                        C(tf) = ptest(tf);
                        F(tf) = ptest(tf);
                        if active(1)== 1
                            Y(tf) = ptest(tf);
                        end
                    end
                end
                if bpfilter == 1
                    fprintf(1,'Masking out areas with patch filter in primary channel...\n');
                    filtersize = side + 2*(bringincr);
                    if finddot == 1 % only filter primary channel. Assume this was channel that was tested to arrive at the filter intensity.
                        O = ordfilt2(Y, 1, ones(filtersize));
                        tf = O>= patchpixel;
                        Y(tf) = ptest(tf);%only replace with 0 if tf is true
                    elseif finddot == 2
                        O = ordfilt2(F, 1, ones(filtersize));
                        tf = O>= patchpixel;
                        F(tf) = ptest(tf);%only replace with 0 if tf is true
                    elseif finddot == 3
                        O = ordfilt2(C, 1, ones(filtersize));
                        tf = O>= patchpixel;
                        C(tf) = ptest(tf);%only replace with 0 if tf is true
                    end
                end
                fprintf(1,'End of all masking routines.\n');
                RY = double( Y );
                RF = double( F );
                RC = double( C );
                point( 5 ) = 5;
                %
                %       i y  j->    Width
                %       | |  x->
                %       v v
                %      Height
                %  avoid the edge
                lowh = frameedge;
                hih  = Height0 - frameedge;
                lowv = frameedge;
                hiv  = Width0  - frameedge;
                sums = zeros( Height0, Width0 );
                initial = zeros( maxpeaks, 3 );
                initialx = zeros( maxpeaks, 3 );
                initialy = zeros( maxpeaks, 3 );
                nearbyZero = zeros(hih,hiv);
                foundspots = 0;
                if( finddot == 1 )
                    fprintf( 1, 'Locating %d peaks,using YFP channel.\n ', maxpeaks );
                elseif( finddot == 2 )
                    fprintf( 1, 'Locating %d peaks, using FRET channel.\n ', maxpeaks );
                else
                    fprintf( 1, 'Locating %d peaks, using CFP channel.\n ', maxpeaks );
                end;
                for kase = finddot
                    if( kase == 1 )
                        M = RY;
                    elseif( kase == 2 )
                        M = RF;
                    else
                        M = RC;
                    end
                    Mcc = zeros(Height0,Width0);
                    Mcr = zeros(Height0,Width0);
                    i = 1:Height0;
                    j = 1:Width0;
                    Mcc(i,j) = cumsum(M,1);
                    Mcr(i,j) = cumsum(Mcc,2);
                    ii = lowh:hih;%avoiding edges
                    jj = lowv:hiv;
                    sums(ii,jj) = Mcr(ii+rindex,jj+rindex) -Mcr(ii+rindex,jj+rindex-side) -(Mcr(ii+rindex-side,jj +rindex) -Mcr(ii+rindex-side,jj+rindex-side));
                    k = 0;
                    got1 = 0;
                    while( k <= maxpeaks && k >= 0)
                        [M1,I] = max(sums);
                        [M2,j] = max(M1);
                        i = I(j);
                        best1y = i; best1x = j; got1 = 1; % get the highest sum in the image
                        best1 = sums(i,j);
                        if best1 > 0.0
                            if best1y-separation < 1
                                best1y = separation + 1;
                            end
                            if best1x-separation <1
                                best1x = separation + 1;
                            end
                            i = (best1y-separation):(best1y+separation);% will zero out this region so not considered during the next pass
                            j = (best1x-separation):(best1x+separation);
                            sums( i, j ) = 0.0;
                            if( k == 0 )
                                k = 1;
                                initial( k, kase ) = best1; initialy( k, kase ) = best1y; initialx( k, kase ) = best1x;
                                foundspots = 1;
                                k = 2;
                            else
                                if( best1 > lowintensity * initial( 1, kase ) )
                                    initial( k, kase ) = best1; initialy( k, kase ) = best1y; initialx( k, kase ) = best1x;
                                    foundspots = foundspots + 1;
                                    k = k + 1;
                                else
                                    k = maxpeaks + 1;
                                end;
                            end;
                        elseif got1 ~= 1
                            foundspots = 0;
                            k = -1;% skip these files
                        else
                            k = -1;
                        end;   % found something
                    end;    % k
                end;      % kase
                fprintf(1, 'Found %d peaks of intensity within the intensity limit. \n', foundspots);
                if foundspots > 0
                    point( 7 ) = 7;
                    rings = zeros( 2*biggestring+1, 1 );
                    for i = 1:2*biggestring+1
                        rings(1, i ) = i;
                    end;
                    AvgRingi = zeros( biggestring+1, maxpeaks, 3 );
                    Totsqr  = zeros( biggestring+1, maxpeaks, 3 );
                    AvgAOIsignal = zeros(maxpeaks, 3);
                    FailedMonoRing = zeros(maxpeaks,3);
                    Monofail = zeros(maxpeaks,3);
                    MinMonToli = zeros(maxpeaks,3);
                    AdjRingRatio = zeros(maxpeaks,3);
                    cxpoint = zeros(maxpeaks,3);
                    cypoint = zeros(maxpeaks,3);
                    FWHM= zeros(3);
                    fitFWHM = zeros(maxpeaks,3);
                    Avgback = zeros(maxpeaks,3);
                    s2n = zeros(maxpeaks,3);
                    for kase = finddot
                        % Make sure peak is the best local maximum. Because of the zeroing out during the
                        % search process, occasionaly the peak should
                        % be moved.
                        if( kase == 1 )
                            M = RY;
                        elseif( kase == 2 )
                            M = RF;
                        else
                            M = RC;
                        end
                        for k = 1:foundspots
                            cxpoint( k, kase ) = initialx( k, finddot );
                            cypoint( k, kase ) = initialy( k, finddot );
                            % if maxmovement is >0 the
                            % local area is defined by maxmovement,
                            % otherwise half the AOI size.
                            if maxmovement == 0
                                wiggle = (side+1)/2;
                            else
                                wiggle = maxmovement;
                            end%
                            altcy = initialy( k, finddot );
                            altcx = initialx( k, finddot );
                            altcyf = altcy;
                            altcxf = altcx;
                            if( side == 1 )
                                altsum = M( altcy, altcx );
                            else
                                altsum = sum(sum(M(altcy-rindex:altcy+rindex,altcx-rindex:altcx+rindex)));
                            end;
                            for i = altcy-wiggle:altcy+wiggle
                                for j = altcx-wiggle:altcx+wiggle
                                    if( side == 1 )
                                        s=M(i,j);
                                    else
                                        s= sum(sum(M(i-rindex:i+rindex,j-rindex:j+rindex)));
                                    end;
                                    if( s > altsum )
                                        altsum = s;
                                        altcyf = i;
                                        altcxf = j;
                                    end;
                                end;
                            end;
                            if( altcx ~= altcxf || altcy ~= altcyf )
                                cxpoint( k, kase ) = altcxf;
                                cypoint( k, kase ) = altcyf;
                                %                                     else
                                %                                         cxpoint( k, kase ) = initialx( k, finddot );
                                %                                         cypoint( k, kase ) = initialy( k, finddot );
                            end;
                            cx = cxpoint( k, kase );
                            cy = cypoint( k, kase );
                            Totsqr( 1, k, kase ) = M( cy, cx );%center pixel of AOI
                            for xring = 2:biggestring+1
                                Totsqr( xring, k, kase ) = 0.0;
                                Totsqr(xring,k,kase) = sum(sum(M(cy-(xring-1):cy+(xring-1),cx-(xring-1):cx+(xring-1))));
                            end;
                            AvgRingi( 1, k, kase ) = Totsqr( 1, k, kase ); %center point
                            for i = 2:biggestring+1
                                AvgRingi( i, k, kase ) = ( Totsqr( i, k, kase ) - Totsqr( i-1, k, kase ) ) / (8*(i-1)); %straight intensity values without background subtraction
                            end
                            AvgAOIsignal(k, kase) = Totsqr(AOIring, k, kase)/side ^ 2;
                        end; %k
                    end;   %kase
                    for kase = 1:3
                        if active(kase) == 1 && kase ~= finddot
                            % Code to move peaks a little bit from image to image because
                            % we are imaging live cells. only performed if
                            % maxmovement >0.
                            if( kase == 1 )
                                M = RY;
                            elseif( kase == 2 )
                                M = RF;
                            else
                                M = RC;
                            end
                            for k = 1:foundspots
                                cxpoint( k, kase ) = cxpoint( k, finddot );
                                cypoint( k, kase ) = cypoint( k, finddot );
                                % If maxmovement is >0 will recenter the peaks in other channels.  Note
                                % this recenters the image outputs.
                                %Then calculates intensity values
                                if maxmovement > 0
                                    altcx = cxpoint( k, finddot );
                                    altcy = cypoint( k, finddot );
                                    altcyf = altcy;
                                    altcxf = altcx;
                                    if( side == 1 )
                                        altsum = M( altcy, altcx );
                                    else
                                        altsum = sum(sum(M(altcy-rindex:altcy+rindex,altcx-rindex:altcx+rindex)));
                                    end;
                                    for i = altcy-maxmovement:altcy+maxmovement
                                        for j = altcx-maxmovement:altcx+maxmovement
                                            if( side == 1 )
                                                s=M(i,j);
                                            else
                                                s= sum(sum(M(i-rindex:i+rindex,j-rindex:j+rindex)));
                                            end;
                                            if( s > altsum )
                                                altsum = s;
                                                altcyf = i;
                                                altcxf = j;
                                            end;
                                        end;
                                    end;
                                    if( altcx ~= altcxf || altcy ~= altcyf )
                                        cxpoint( k, kase ) = altcxf;
                                        cypoint( k, kase ) = altcyf;
                                    end;
                                end
                                cx = cxpoint( k, kase );
                                cy = cypoint( k, kase );
                                Totsqr( 1, k, kase ) = M( cy, cx );%center pixel of AOI
                                for xring = 2:biggestring+1
                                    Totsqr( xring, k, kase ) = 0.0;
                                    Totsqr(xring,k,kase) = sum(sum(M(cy-(xring-1):cy+(xring-1),cx-(xring-1):cx+(xring-1))));
                                end;
                                AvgRingi( 1, k, kase ) = Totsqr( 1, k, kase ); %center point
                                for i = 2:biggestring+1
                                    AvgRingi( i, k, kase ) = ( Totsqr( i, k, kase ) - Totsqr( i-1, k, kase ) ) / (8*(i-1)); %straight intensity values without background subtraction
                                end
                                AvgAOIsignal(k, kase) = Totsqr(AOIring, k, kase)/side ^ 2;
                            end; %k
                        end; %active kases only
                    end;   %kase
                    point( 8 ) = 8;
                    if( debug ==  1 )
                        fprintf( 1, 'got past point( 8 )\n' );
                    end;
                    fprintf(1, 'Testing peaks to see if they meet search criteria.\n');
                    mtol( 1 )   = YFPmonotonicity;  FWHM( 1 ) = YFPFWHM;  s2ntarget( 1 ) = YFPsignal2noise;
                    mtol( 2 )   = FRETmonotonicity; FWHM( 2 ) = FRETFWHM; s2ntarget( 2 ) = FRETsignal2noise;
                    mtol( 3 )   = CFPmonotonicity;  FWHM( 3 ) = CFPFWHM;  s2ntarget( 3 ) = CFPsignal2noise;
                    point( 9 ) = 9;
                    if( debug == 1 )
                        fprintf( 1, 'got past point( 9 )\n' );
                    end;
                    dot = 1;
                    while( dot <= foundspots )
                        keep = 1;
                        for kase = 1:3
                            if active(kase)== 1
                                if kase == 1
                                    channel = 'YFP';
                                elseif kase == 2
                                    channel = 'FRET';
                                else
                                    channel = 'CFP';
                                end
                                Backpixs = (( (2 * lastbring) + 1 ) ^ 2 )- (( (2 * firstbring) - 1 ) ^ 2);
                                Avgback( dot, kase ) = ( Totsqr(lastbring+1,dot,kase) - Totsqr(firstbring,dot,kase) ) / Backpixs;
                                s2n(dot,kase) = AvgAOIsignal(dot,kase) / Avgback(dot,kase);
                                if kase == finddot || (kase ~= finddot && limit1 == 0)
                                    if( (( AvgAOIsignal( dot, kase ) - Avgback( dot, kase ) ) < 0.0) && ( active( kase ) == 1 ) )
                                        keep = failaction;
                                        fprintf(    1, '      FAILED in %s channel. Signal(%6.2f) < background(%6.2f).\n',channel, AvgAOIsignal(dot,kase),Avgback(dot,kase));
                                    end;
                                end
                            end
                            if kase == finddot || (kase ~= finddot && limit1 == 0)
                                if( active( kase ) == 1 && keep == 1 )         % only active channels are tested
                                    if AvgAOIsignal(dot,kase)==0 || Avgback(dot,kase)== 0
                                        keep = failaction;
                                        fprintf(    1, '     FAILED in %s channel. A too bright area was set to 0.\n', channel );
                                    end;
                                end;
                            end
                            if kase == finddot || (kase ~= finddot && limit1 == 0)
                                if( active( kase ) == 1 && keep == 1 )         % only active channels are tested
                                    if AvgAOIsignal(dot,kase) /Avgback(dot,kase) < s2ntarget(kase)
                                        keep = failaction;
                                        fprintf(    1, '     FAILED in %s channel.  Signal to noise too low. Target = %5.2f, Observed = %5.2f.\n', channel, s2ntarget(kase), s2n(dot,kase) );
                                    end;
                                end;
                            end
                            % monotonic
                            if( active( kase ) == 1 && keep == 1 )
                                monotonic = 1;
                                for i = 1:AOIring
                                    AdjRingRatio(dot,kase,i) = AvgRingi( i+1, dot, kase )/AvgRingi( i, dot, kase )-1;
                                    MinMonToli(dot,kase,i)= AdjRingRatio(dot,kase,i);
                                    if kase == finddot || (kase ~= finddot && limit1 == 0)
                                        if ((1 + mtol(kase)) * AvgRingi( i, dot, kase ) <  AvgRingi( i+1, dot, kase )) && monotonic == 1 %mtol is monotonicity tolerance
                                            monotonic = 0;
                                            FailedMonoRing(dot,kase) = i;
                                            Monofail(dot,kase) = ceil(((AvgRingi( i+1, dot, kase )/AvgRingi( i, dot, kase ))-1)*100)/100;%would require this mtol to pass
                                        end
                                    end
                                end
                                if( monotonic == 0 )
                                    keep = failaction;
                                    fprintf(    1, '     FAILED in the %s channel. Failed at ring %d. Target monotonicity tolerance = %5.2f. Would require %5.2f to pass.\n', channel,FailedMonoRing(dot,kase),mtol( kase ),Monofail(dot,kase)  );
                                end;
                            end
                            % active and keep
                            %  Next bit of code checks whether peak is gaussian and meets the input widths
                            if( active( kase ) == 1 && keep == 1 );
                                start( 1 ) = FWHM(kase);
                                posit = zeros(1,(2*biggestring + 1));
                                for i = 1:biggestring +1
                                    z = i + biggestring;
                                    posit( z ) = AvgRingi( i,dot, kase )-Avgback(dot,kase);
                                end;
                                for i = 1:biggestring
                                    posit(i)= posit((2*biggestring)+2-i);
                                end
                                FWHMconst = 2*sqrt(2*log(2));
                                options = optimset('MaxFunEvals',300,'MaxIter',300);
                                [c,fval,exitflag,output] = fminsearch( @gaussian, start,options );
                                if exitflag(1)< 1
                                   if kase == finddot || (kase ~= finddot && limit1 == 0)
                                    keep = failaction;
                                    fprintf(    1, '     FAILED in %s channel during Gaussian fit. Did not converge.\n', channel );
                                   else
                                     fprintf(    1, '    During Gaussian fit %s channel did not converge. AOI kept because limits off for this channel.\n', channel );   
                                   end
                                end
                                if exitflag(1) == 0
                                    fprintf(    1, '     Maximum number of function evaluations or iterations was reached.\n     Function evals %d. Iterations %d\n', output(2),output(3));
                                end
                                if exitflag(1) == 1
                                    fitFWHM(dot,kase) = abs(c(1)* FWHMconst);
                                    if kase == finddot || (kase ~= finddot && limit1 == 0)
                                        if( fitFWHM( dot,kase ) > FWHM( kase ) )  % peak is too broad
                                            keep = failaction;
                                            fprintf(    1, '     FAILED in %s channel during Gaussian fit. Target FWHM= %5.2f. Observed FWHM= %5.2f. \n', channel, FWHM( kase ),fitFWHM( dot,kase ) );
                                        end;
                                        if( fitFWHM( dot,kase ) < FWHMminimum )  % too fast
                                            keep = failaction;
                                            fprintf(    1, '     FAILED in %s channel. Default FWHMminimum = %5.2f. Observed FWHM= %5.2f. Probably hot pixel.\n', channel, FWHMminimum, fitFWHM( dot,kase ) );
                                        end;
                                    end
                                end
                            end;  % active and keep
                        end;   %kase
                        if( keep == 1 )
                            MinMonTol = max(MinMonToli,[],3);
                            if( writesmalltiff == 1 )
                                Pad2 = context; %Pad the image again to accomodate the context in the output tiffs
                                subwidth  = 2 * (Pad2 + biggestring);
                                subheight = 2 * (Pad2 + biggestring);
                                if active(1) == 1
                                    padY = padarray(Y, [ Pad2 Pad2 ], MinY);
                                    startx = cxpoint( dot, 1 ) - biggestring;
                                    starty = cypoint( dot, 1 ) - biggestring;
                                    suby = imcrop( padY, [ startx starty subwidth subheight ] );  % order of acquisition
                                    imwrite( suby, strcat( datadir, dirsep, num2str(dot,'%03.3d'),Q{1,e}), 'tif','Compression','none' );
                                end
                                if active(2) == 1
                                    startx = cxpoint( dot, 2 ) - biggestring;
                                    starty = cypoint( dot, 2 ) - biggestring;
                                    padF = padarray(F, [ Pad2 Pad2], MinF);
                                    subf = imcrop( padF, [ startx starty subwidth subheight ] );
                                    imwrite( subf, strcat( datadir, dirsep, num2str(dot,'%03.3d'),Q{2,e}), 'tif','Compression','none' );
                                end
                                if active(3) == 1
                                    startx = cxpoint( dot, 3 ) - biggestring;
                                    starty = cypoint( dot, 3 ) - biggestring;
                                    padC = padarray(C, [ Pad2 Pad2 ], MinC);
                                    subc = imcrop( padC, [ startx starty subwidth subheight ] );
                                    imwrite( subc, strcat( datadir, dirsep, num2str(dot,'%03.3d'),Q{3,e}), 'tif', 'Compression','none' );
                                end
                                if active(4) == 1
                                    startx = cxpoint( dot, finddot ) - biggestring;
                                    starty = cypoint( dot, finddot ) - biggestring;
                                    padD = padarray(D, [ Pad2 Pad2 ], min(min(D)));
                                    subd = imcrop( padD, [ startx starty subwidth subheight ] );
                                    imwrite( subd, strcat( datadir, dirsep, num2str(dot,'%03.3d'),Q{4,e}), 'tif', 'Compression','none' );
                                end
                            end;
                            filename = strcat(num2str(dot,'%03.3d'),Q{5,e});
                            fprintf( 1, 'WRITTEN to %s\n', filename);
                            adjustYFP = AvgAOIsignal(dot, 1) - Avgback(dot,1);
                            adjustFRET = AvgAOIsignal(dot,2) - Avgback(dot,2);
                            adjustCFP = AvgAOIsignal(dot,3) - Avgback(dot,3);
                            Spillover_tot  = YFPspill * adjustYFP + CFPspill * adjustCFP;%if calculating spillover factors this just returns the adjusted CFP or YFP
                            fret       = adjustFRET / Spillover_tot; % if calculating spillover factors this variable represents the spillover factor
                            fprintf( fjmp, '%-20s %10s %10s %10s %2d %8.3f %8.3f  %8.3f %8.3f  %8.3f %8.3f %8.3f %8.3f %8.3f  %8.3f %8.3f',...
                                filename,strain,donor,acceptor,dot, ...
                                AvgAOIsignal(dot,1),Avgback(dot,1),AvgAOIsignal(dot,2),Avgback(dot,2),AvgAOIsignal(dot,3),Avgback(dot,3), ...
                                adjustYFP, adjustFRET, adjustCFP, Spillover_tot, fret);
                            for kase = 1:3
                                if( active( kase ) == 1 )
                                    fprintf( fjmp, '  %8.3f  %8.3f  %8.3f  %d  %d ', fitFWHM(dot,kase), s2n(dot,kase),MinMonTol(dot,kase),cxpoint(dot,kase),cypoint(dot,kase) );
                                else
                                    fprintf( fjmp, '  %8.3f  %8.3f  %8.3f %d  %d ', 0.0, 0.0, 0.0, 0, 0 );
                                end;
                            end;
                            fprintf( fjmp, '\n' );
                        end;
                        dot = dot + 1;
                    end; % dot
                    point( 10 ) = 10;
                    fail = 0;
                end;
            end%skipall
        end
        toc
    end
    fclose( fjmp );
end
beep

function[lsr,gfit] = gaussian(c)
%least squares fit to a Gaussian distribution
global biggestring posit

mid = biggestring +1;
range = 2*biggestring +1;
ht = posit(mid);
gfit = zeros(1,range);
for i = 1:range
    z = i-mid;%mean = 0
    gfit(i) = ht*(1/(exp(((z)^2)/(2*(c^2)))));
end
lsr = 0;
for i = 1:range
    lsr = lsr + ( gfit(i) - (posit( i )) )^2;
end

Contact us