home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / FINDMAX2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  112 lines
  1. function [gt,d]=findmax2(data, flag)
  2. %FINDMAX2 Interpolates the maxima in a matrix of data.
  3. %
  4. %    Function used for returning all the maxima in a set of
  5. %    data (matrices). The maxima are calculated by
  6. %    quadratic interpolation.
  7.  
  8. if nargin < 2, flag = 0, end
  9. % 2-D
  10. d = [1.5, 1.5];
  11. [m,n]=size(data);
  12. % Factor for error control:    
  13. factor = 5;  % Keep error five times less than nearest constraint
  14.  
  15. % Constants
  16. diagix = [-1 -1 1 1];
  17. diagix2 = [1 -1 -1 1];
  18. % Use point [0,0], [1,0], [0,1], [0, -1], [-1 0] and [1,1] on grid:
  19. diagix3 = [0 1 0 0 -1 1 ];
  20. diagix4 = [0 0 1 -1 0 1 ];
  21. % Point [-1, -1] used to check goodness of quadratic fit. 
  22. xerr = [-1; -1];
  23.  
  24. invH = [  ... 
  25. -1.0000    0.5000         0         0    0.5000         0
  26. 0.5000   -0.5000   -0.5000         0         0    0.5000
  27. -1.0000         0    0.5000    0.5000         0         0
  28. 0    0.5000         0         0   -0.5000         0
  29. 0         0    0.5000   -0.5000         0         0
  30. 1.0000         0         0         0         0         0
  31. ];
  32.  
  33. % Find the peaks
  34. crmax = (data>[data(1,:)-1;data(1:m-1,:)]) & (data>=[data(2:m,:);data(m,:)-1]) & ...
  35.     (data>[data(:,1)-1,data(:,1:n-1)]) & (data>=[data(:, 2:n),data(:,n)-1]);
  36.  
  37. f = find(crmax);
  38. for i = f'
  39.     ni = 1 + floor((i-0.5)/m);
  40.     mi = i - (ni-1)*m;
  41.     ix = mi - diagix;
  42.     iy = ni - diagix2;
  43.     fix = find(ix > 0 & ix <= m & iy > 0 & iy <= n); 
  44.     if any(data(mi, ni) <= data(ix(fix), iy(fix)))
  45.         crmax(mi, ni) = 0; 
  46.     end
  47. end
  48. f = find(crmax)';
  49. H =zeros(2,2); 
  50. % Fit a quadratic
  51. cnt = 0; 
  52. for i = f
  53.     cnt = cnt + 1;
  54.     ni = 1 + floor((i-0.5)/m);
  55.     mi = i - (ni-1)*m;
  56.     if (mi > 1 & mi < m & ni > 1 & ni < n) 
  57.         ix = mi - diagix3;
  58.         iy = ni - diagix4;
  59.         for k = 1:6
  60.             y(k,1) = data(ix(k),iy(k)); 
  61.         end
  62.         coeffs = invH*y; 
  63.         H(1,1) = coeffs(1);
  64.         H(1,2) = coeffs(2);
  65.         H(2,1) = H(1,2);
  66.         H(2,2) = coeffs(3); 
  67.         b = coeffs(4:5);
  68.         c = coeffs(6);
  69.         x = -0.5 * (H\b); 
  70.         gt(cnt) = x'*H*x + b'*x + c;
  71. % Error control 
  72.         if flag 
  73.             err2 = abs(xerr'*H*xerr + b'*xerr + c - data(mi+1, ni+1 ));
  74.             mult = 1; 
  75.             if err2 < 1/20*abs(gt(cnt))
  76.                 mult = 1.25;
  77.             elseif err2 > (5*abs(gt(cnt))+1e-5)
  78.                 mult = 0.75;
  79.             end
  80.  
  81.             pky = quad2(data(mi-1:mi+1,ni));
  82.             f = abs((pky + 1.1e-5) /(factor *(data(mi,ni)-pky) + 1e-5));
  83.             d(1,2) = min(d(1,2), mult * ((f>=0.3) + (0.7 + f)*(f< 0.3) +(f>1)*(log(f)/100)));
  84.             pkx =  quad2(data(mi, ni-1:ni+1)');
  85.             f = abs((pkx + 1.1e-5) /(factor  *(data(mi,ni) - pkx) + 1e-5));
  86.             d(1,1) = min(d(1,1), mult * ((f>=0.3) + (0.7 + f)*(f< 0.3) +(f>1)*(log(f)/100)));
  87.         end
  88.     elseif mi > 1 & mi < m
  89.         gt(cnt) = quad2(data(mi-1:mi+1,ni));
  90.         if flag 
  91.             pky = gt(cnt);
  92.             f = abs((pky + 1.1e-5) /(factor *(data(mi,ni) - pky) + 1e-5));
  93.             d(1,2) = min(d(1,2), (f>=0.3) + (0.7 + f)*(f< 0.3) +(f>1)*(log(f)/100));
  94.             if (pky > -1e-5), d(1,1) = min(1, d(1,1)); end
  95.         end
  96.     elseif ni > 1 & ni < n
  97.         gt(cnt) = quad2(data(mi, ni-1:ni+1)');
  98.         if flag 
  99.             pkx = gt(cnt); 
  100.             f = abs((pkx + 1.1e-5) /(factor  *(data(mi,ni) - pkx) + 1e-5));
  101.             d(1,1) = min(d(1,1), (f>=0.3) + (0.7 + f)*(f< 0.3) +(f>1)*(log(f)/100));
  102.             if (pkx > -1e-5), d(1,2) = min(1, d(1,2)); end
  103.         end
  104.     else
  105.         gt(cnt) = data(mi, ni);
  106.         if (gt(cnt) > -1e-5), 
  107.             d(1,1) = min(1, d(1,1)); 
  108.             d(1,2) = min(1, d(1,2)); 
  109.         end
  110.     end
  111. end
  112.