How can I find the maximum gradient around a data point?

I want to find the x/y indices on a surface where the gradient reaches its maximum around a point P.

For example, let's take the following surface Z:

[X,Y] = meshgrid(-2:.2:2);
Z = X.*exp(-X.^2 - Y.^2);
Z(Z<0) = 0;
[DX,DY] = gradient(Z,.2,.2);
g = sqrt(DX.^2+DY.^2);

figure
contour(X,Y,Z)
hold on
quiver(X,Y,DX,DY)
hold off

Let's say the position of the peak, Z(21,28), is known. Now, I'm looking for the x/y indices, where g (gradient) reaches its highest values around this peak, which would be an oval-like line. How can it be done?

Answer 1

[Soheil Solhjoo]

I tried the following, which does not produce the best results, but works along the idea.

    %%Scan g from the selected point P(y,x), and detect the border
    % border : where g reaches its maximum somewhere close to P
    P = [21,28]; % selected point: the peak, in this case
    g_x = g(P(1,1),:);
    r_b = right_border_x(g_x,P(1,2)); % see below
    r0 = r_b - P(1,2); % initial guess for the radius
    range = 0; % this is the range for function track_border_r (see below)
    tracked_border = track_border_r(grad_temp,P(1,:),r0,360,range);

For running the code, I wrote two functions. "x_right_border":

    function x_right_border = right_border_x(gradient_vector,position_x)
    % The function will search from position_x to the right hand side, i.e.
    % values larger than position_x, till reach the border, defined via the
    % maximum gradient.
    g_xo = gradient_vector(position_x); % Gradient at the original position_x
    g_xr = gradient_vector(position_x+1); % Gradient at its right neighbor
    if g_xr < g_xo
        x_right_border = position_x;
    else
        position_x = position_x+1;
        x_right_border = right_border_x(gradient_vector,position_x);
    end
    end

, and "track_border_r"

    function r_tracked = track_border_r(gradient_vector,center,r0,thetha_inc,range)
    % gradient_vector : gradient of the surcface
    % center(y,x) : selected center
    % r0 : initial value of r for theta=0
    % thetha_inc : number of increments of thetha, between 0 and 2*pi
    % range: range for searching for max in the gradient vector around center
    x_c = center(1,2);
    y_c = center(1,1);
    i = 1;
    for theta = 0 : 2*pi/thetha_inc : 2*pi*(1-1/thetha_inc)
        [x_i,y_i] = pol2cart(theta,r0);
        x_r = x_i + x_c;
        x_r = round(x_r-range:x_r+range);
        y_r = round(y_i + y_c);
        y_r = round(y_r-range:y_r+range);
        g_temp = gradient_vector(y_r,x_r);
        [~,ind] = max(g_temp(:));
        [ind_row, ind_col] = ind2sub(size(g_temp),ind);
        x_max = x_r(1) + ind_col - 1;
        y_max = y_r(1) + ind_row - 1;
        r0 = sqrt((x_max-x_c)^2 + (y_max-y_c)^2);
        tracked(i,:) = [y_max,x_max];
        i = i + 1;
    end
    r_tracked = unique(tracked,'rows');
    end

By changing the values of "range" and "thetha_inc" in function "track_border_r" its accuracy can be controlled to some extent; however, the code does not work flawlessly.