¿Alguien sabe cómo hacer un contorno activo abierto? Estoy familiarizado con los contornos cerrados y tengo varios programas de Matlab que los describen. Intenté cambiar la matriz P pero eso no fue suficiente: mi código de Matlab es el siguiente:

% Read in the test image
GrayLevelClump=imread('cortex.png'); 
cortex_grad=imread('cortex_grad.png');
rect=[120 32 340 340];
img=imcrop(GrayLevelClump,rect);
grad_mag=imcrop(cortex_grad(:,:,3),rect);
% Draw the initial snake as a line
a=300;b=21;c=21;d=307;
snake = brlinexya(a,b,c,d);    % startpoints,endpoints
x0=snake(:,1);
y0=snake(:,2);
N=length(x0);
x0=x0(1:N-1)';
y0=y0(1:N-1)';

% parameters for the active contour
alpha = 5;
beta = 10;
gamma = 1;
kappa=1;
iterations = 50;
% Make a force field
[u,v] = GVF(-grad_mag, 0.2, 80);
 %disp(' Normalizing the GVF external force ...');
mag = sqrt(u.*u+v.*v);
px = u./(mag+(1e-10)); py = v./(mag+(1e-10));
%  Make the coefficient matrix P
x=x0;
y=y0;
N = length(x0);
a = gamma*(2*alpha+6*beta)+1;
b = gamma*(-alpha-4*beta);
c = gamma*beta;
P = diag(repmat(a,1,N));
P = P + diag(repmat(b,1,N-1), 1) + diag(   b, -N+1);
P = P + diag(repmat(b,1,N-1),-1) + diag(   b,  N-1);
P = P + diag(repmat(c,1,N-2), 2) + diag([c,c],-N+2);
P = P + diag(repmat(c,1,N-2),-2) + diag([c,c], N-2);
P = inv(P);

d = gamma * (-alpha);
e = gamma * (2*alpha);
% Do the modifications to the matrix in order to handle OAC
P(1:2,:) = 0;
P(1,1) = -gamma;
P(2,1) = d;
P(2,2) = e;
P(2,3) = d;

P(N-1:N,:) = 0;
P(N-1,N-2) = d;
P(N-1,N-1) = e;
P(N-1,N) = d;
P(N,N) = -gamma;

x=x';
y=y'; 
figure,imshow(grad_mag,[]);
hold on, plot([x;x(1)],[y;y(1)],'r');
plot(x([1 end]),y(([1 end])), 'b.','MarkerSize', 20);

for ii = 1:50
   % Calculate external force
   vfx = interp2(px,x,y,'*linear');
   vfy = interp2(py,x,y,'*linear');
   tf=isnan(vfx);
    ind=find(tf==1);
    if ~isempty(ind)
        vfx(ind)=0;
    end
    tf=isnan(vfy);
    ind=find(tf==1);
    if ~isempty(ind)
        vfy(ind)=0;
    end

   % Move control points
   x = P*(x+gamma*vfx);
   y = P*(y+gamma*vfy);
   %x = P * (gamma* x + gamma*vfx);
   %y = P * (gamma* y + gamma*vfy);
    if mod(ii,5)==0
      plot([x;x(1)],[y;y(1)],'b')
    end
end
plot([x;x(1)],[y;y(1)],'r')

Modifiqué la matriz de coeficientes P para manejar casos con contornos activos abiertos, pero eso claramente no es suficiente.

Mi imagen: