/normxcorr/trunk

To get this branch, use:
bzr branch http://suren.me/webbzr/normxcorr/trunk

« back to all changes in this revision

Viewing changes to findpeak.m

  • Committer: Suren A. Chilingaryan
  • Date: 2009-01-15 13:50:29 UTC
  • Revision ID: csa@dside.dyndns.org-20090115135029-wleapicg9a4593tp
Initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
findpeak.m
 
1
function [xpeak, ypeak, max_f] = findpeak(f,subpixel)
 
2
%FINDPEAK Find extremum of matrix.
 
3
%   [XPEAK,YPEAK,MAX_F] = FINDPEAK(F,SUBPIXEL) finds the extremum of F,
 
4
%   MAX_F, and its location (XPEAK, YPEAK). F is a matrix. MAX_F is the maximum
 
5
%   absolute value of F, or an estimate of the extremum if a subpixel
 
6
%   extremum is requested.
 
7
%
 
8
%   SUBPIXEL is a boolean that controls if FINDPEAK attempts to estimate the
 
9
%   extremum location to subpixel precision. If SUBPIXEL is false, FINDPEAK
 
10
%   returns the coordinates of the maximum absolute value of F and MAX_F is
 
11
%   max(abs(F(:))). If SUBPIXEL is true, FINDPEAK fits a 2nd order
 
12
%   polynomial to the 9 points surrounding the maximum absolute value of
 
13
%   F. In this case, MAX_F is the absolute value of the polynomial evaluated
 
14
%   at its extremum.
 
15
%
 
16
%   Note: Even if SUBPIXEL is true, there are some cases that result
 
17
%   in FINDPEAK returning the coordinates of the maximum absolute value
 
18
%   of F:
 
19
%   * When the maximum absolute value of F is on the edge of matrix F.
 
20
%   * When the coordinates of the estimated polynomial extremum would fall
 
21
%     outside the coordinates of the points used to constrain the estimate.
 
22
 
 
23
%   Copyright 1993-2004 The MathWorks, Inc.
 
24
%   $Revision $  $Date: 2004/10/20 17:54:47 $
 
25
 
 
26
% get absolute peak pixel
 
27
[max_f, imax] = max(abs(f(:)));
 
28
[ypeak, xpeak] = ind2sub(size(f),imax(1));
 
29
    
 
30
if ~subpixel || ...
 
31
    xpeak==1 || xpeak==size(f,2) || ypeak==1 || ypeak==size(f,1) % on edge
 
32
    return % return absolute peak
 
33
    
 
34
else
 
35
    % fit a 2nd order polynomial to 9 points  
 
36
    % using 9 pixels centered on irow,jcol    
 
37
    u = f(ypeak-1:ypeak+1, xpeak-1:xpeak+1);
 
38
    u = u(:);
 
39
    x = [-1 -1 -1  0  0  0  1  1  1]';
 
40
    y = [-1  0  1 -1  0  1 -1  0  1]';    
 
41
 
 
42
    % u(x,y) = A(1) + A(2)*x + A(3)*y + A(4)*x*y + A(5)*x^2 + A(6)*y^2
 
43
    X = [ones(9,1),  x,  y,  x.*y,  x.^2,  y.^2];
 
44
    
 
45
    % u = X*A
 
46
    A = X\u;
 
47
 
 
48
    % get absolute maximum, where du/dx = du/dy = 0
 
49
    x_offset = (-A(3)*A(4)+2*A(6)*A(2)) / (A(4)^2-4*A(5)*A(6));
 
50
    y_offset = -1 / ( A(4)^2-4*A(5)*A(6))*(A(4)*A(2)-2*A(5)*A(3));
 
51
 
 
52
    if abs(x_offset)>1 || abs(y_offset)>1
 
53
        % adjusted peak falls outside set of 9 points fit,
 
54
        return % return absolute peak
 
55
    end
 
56
    
 
57
    % return only one-tenth of a pixel precision
 
58
    x_offset = round(10*x_offset)/10;
 
59
    y_offset = round(10*y_offset)/10;    
 
60
    
 
61
    xpeak = xpeak + x_offset;
 
62
    ypeak = ypeak + y_offset;    
 
63
    
 
64
    % Calculate extremum of fitted function
 
65
    max_f = [1 x_offset y_offset x_offset*y_offset x_offset^2 y_offset^2] * A;
 
66
    max_f = abs(max_f);
 
67
    
 
68
end