off-motion-receptive-fields/receptiveFields/fitDifferenceOfGaussiansFromNeuronStimTables.m

function flyFits = fitDifferenceOfGaussiansFromNeuronStimTables(StimTableNeuronStim)
tuningMethod = 'extreme';
flyFits = {};
for neuronTypeInd = 1: size(StimTableNeuronStim, 1)
    for stimInd = 2: size(StimTableNeuronStim, 2)
        TempTable = StimTableNeuronStim{neuronTypeInd, stimInd};
        for flyInd = 1: size(TempTable, 1)
            tuningInd = find(strcmp(TempTable.nonParamTuning{1}.tuningMethods, tuningMethod));
            tuningCurve = squeeze(TempTable.nonParamTuning{flyInd}.tc{tuningInd});
            if stimInd == 3 || stimInd == 5
                tuningCurve = tuningCurve(:, 3: end - 2);
            end
            if stimInd > 3
                gaussianSign = -1;
            else
                gaussianSign = 1;
            end
            % Still need to crop extremes for consistency with normal
            % fitting.
            xArray = 1: size(tuningCurve, 2);
            fitStruct(size(tuningCurve, 1), 1).fit = [];
            fitStruct(size(tuningCurve, 1), 1).gof = [];
            fitStruct(size(tuningCurve, 1), 1).info = [];
            for kCell = 1: size(tuningCurve, 1)
                [fitStruct(kCell).fit, ...
                 fitStruct(kCell).gof.rsquare] = fitDifferenceOfGaussianConstraint(...
                                                    xArray, ...
                                                    tuningCurve(kCell, :), ...
                                                    gaussianSign);
                figure;
                scatter(xArray, tuningCurve(kCell, :), 'filled');
                hold on;
                plot(xArray, differenceOfGaussians1D(fitStruct(kCell).fit, xArray));
                pause(0.01)
            end
            flyFits{neuronTypeInd, stimInd - 1, flyInd} = fitStruct;
            clear fitStruct;
            close all;
        end
    end
end

end