fitDifferenceOfGaussianConstraint.m 1.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647
  1. function [a, rSquared] = fitDifferenceOfGaussianConstraint(x, y, gaussianSign)
  2. %%
  3. % DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
  4. % a(4) * exp(-((x - a(2)) / a(5)).^2);
  5. % lsq = @(a) sum(F(a) - y).^2;
  6. DoG = @(a) differenceOfGaussians1D(a, x);
  7. % DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
  8. % a(4) * exp(-((x - a(2)) / a(5)).^2);
  9. lsq = @(a) sum((DoG(a) - y).^2);
  10. [minY, argMinY] = min(y);
  11. [maxY, argMaxY] = max(y);
  12. xStep = mean(diff(x));
  13. positiveAmpBounds = [0, 2 * range(y)];
  14. negativeAmpBounds = [- 2 * range(y), 0];
  15. positionBounds = x([1, end]);
  16. if gaussianSign > 0
  17. centerAmpBounds = positiveAmpBounds;
  18. centerWidthBounds = [xStep, range(x)];
  19. surroundAmpBounds = negativeAmpBounds;
  20. surroundWidthBounds = [2 * xStep, range(x)];
  21. a0 = [maxY, argMaxY, xStep, minY, range(x) / 2];
  22. else
  23. centerAmpBounds = negativeAmpBounds;
  24. centerWidthBounds = [2 * xStep, range(x)];
  25. surroundAmpBounds = positiveAmpBounds;
  26. surroundWidthBounds = [2 * xStep, range(x)];
  27. a0 = [minY, argMinY, xStep, maxY, range(x) / 2];
  28. end
  29. % center should be narrower than surround: a(3) < a(5) <=> a(3) - a(5) < 0
  30. A = [0 0 3 0 -5];
  31. b = 0;
  32. Aeq = [];
  33. beq = [];
  34. lb = [centerAmpBounds(1), positionBounds(1), centerWidthBounds(1), ...
  35. surroundAmpBounds(1), surroundWidthBounds(1)];
  36. ub = [centerAmpBounds(2), positionBounds(2), centerWidthBounds(2), ...
  37. surroundAmpBounds(2), surroundWidthBounds(2)];
  38. [a, fitResidue] = fmincon(lsq, a0, A, b, Aeq, beq, lb, ub);
  39. rSquared = 1 - fitResidue / sum((y - mean(y)).^2);
  40. end