Code covered by the BSD License  

Highlights from
FRETSCAL

image thumbnail

FRETSCAL

by

 

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