1234567891011121314151617181920212223242526272829303132333435363738394041424344454647 |
- function [a, rSquared] = fitDifferenceOfGaussianConstraint(x, y, gaussianSign)
- %%
- % DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
- % a(4) * exp(-((x - a(2)) / a(5)).^2);
- % lsq = @(a) sum(F(a) - y).^2;
- DoG = @(a) differenceOfGaussians1D(a, x);
- % DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
- % a(4) * exp(-((x - a(2)) / a(5)).^2);
- lsq = @(a) sum((DoG(a) - y).^2);
- [minY, argMinY] = min(y);
- [maxY, argMaxY] = max(y);
- xStep = mean(diff(x));
- positiveAmpBounds = [0, 2 * range(y)];
- negativeAmpBounds = [- 2 * range(y), 0];
- positionBounds = x([1, end]);
- if gaussianSign > 0
- centerAmpBounds = positiveAmpBounds;
- centerWidthBounds = [xStep, range(x)];
- surroundAmpBounds = negativeAmpBounds;
- surroundWidthBounds = [2 * xStep, range(x)];
- a0 = [maxY, argMaxY, xStep, minY, range(x) / 2];
- else
- centerAmpBounds = negativeAmpBounds;
- centerWidthBounds = [2 * xStep, range(x)];
- surroundAmpBounds = positiveAmpBounds;
- surroundWidthBounds = [2 * xStep, range(x)];
- a0 = [minY, argMinY, xStep, maxY, range(x) / 2];
- end
- % center should be narrower than surround: a(3) < a(5) <=> a(3) - a(5) < 0
- A = [0 0 3 0 -5];
- b = 0;
- Aeq = [];
- beq = [];
- lb = [centerAmpBounds(1), positionBounds(1), centerWidthBounds(1), ...
- surroundAmpBounds(1), surroundWidthBounds(1)];
- ub = [centerAmpBounds(2), positionBounds(2), centerWidthBounds(2), ...
- surroundAmpBounds(2), surroundWidthBounds(2)];
- [a, fitResidue] = fmincon(lsq, a0, A, b, Aeq, beq, lb, ub);
- rSquared = 1 - fitResidue / sum((y - mean(y)).^2);
- end
|