No BSD License  

Select ROI in image using spline

by

 

04 Oct 2006 (Updated )

A usefull program to select smooth region in image

mask=roispline(I,kind,tension)
function mask=roispline(I,kind,tension)

% function [mask,perimeter,area]=roispline(I,str,tension)
%
% INPUT :
% - I       : Grayscale or color image
% - kind    : String specifying the kind of spline ('natural' for natural 
%             cubic spline, 'cardinal' for cardinal cubic spline). 
%             DEFAULT = 'natural'
% - tension : Tension parameter between 0 and 1 for cardinal spline.
%             DEFAULT = 0
%
% OUTPUT :
% - mask    : Binary mask of segmented ROI


warning off

if nargin<2
    kind='natural';
elseif nargin<3
    tension=0;
end

siz=size(I);
imshow(I);
title('Click right button or close curve clicking near first point');
hold on;


% initializes number of points for each step
npoints=sqrt(siz(1)*siz(2))*.1;
dim=10/835*sqrt(siz(1)*siz(2));
if dim<8
    dim=8;
end

bottone=1;
i=0;
flag=1;
points=[];
switch kind
    case 'natural'
        while flag                   
            i=i+1;
            [cor(1,i),cor(2,i),bottone]=ginput(1);
            if bottone ~= 1
                if i>1
                    cor(:,i)=cor(:,1);
                    flag=0;
                else
                    break;
                end
            end
            plot(cor(1,i),cor(2,i),'y.');
    
            if i>1      
                if (norm(cor(:,1)-cor(:,i))<dim & i>2) & flag
                    cor(:,i)=cor(:,1);
                    flag=0;
                end
                imshow(I)
                hold on;
                spcv = cscvn(cor);              
                points=myfnplt(spcv,npoints*i);
                plot(points(1,:),points(2,:),'b','LineWidth',2);
                plot(cor(1,:),cor(2,:),'y.');
                drawnow;
            end
        end
    case 'cardinal'                  
        while flag                         
            i=i+1;
            [x(i),y(i),bottone]=ginput(1);
            if bottone ~= 1
                if i>1
                    x(i)=[];
                    y(i)=[];
                    Px=[x(end) x x(1) x(2)];
                    Py=[y(end) y y(1) y(2)];
                    i=i-1;
                    flag=0;
                else
                    break;
                end
            end            
            plot(x(i),y(i),'y.');  
            
            if i>1 
                if flag
                    if (norm([x(1) y(1)]-[x(i) y(i)])<dim & i>2)
                        x(i)=[];
                        y(i)=[];
                        Px=[x(end) x x(1) x(2)];
                        Py=[y(end) y  y(1) y(2)];
                        flag=0;
                    else
                        Px=[x(1) x x(end)];
                        Py=[y(1) y y(end)];
                    end
                end
                pointsx=[];
                pointsy=[];
                for k=1:length(Px)-3
                    [xvec,yvec]=EvaluateCardinal2DAtNplusOneValues([Px(k),Py(k)],[Px(k+1),Py(k+1)],[Px(k+2),Py(k+2)],[Px(k+3),Py(k+3)],tension,npoints);
                    pointsx=[pointsx, xvec];
                    pointsy=[pointsy, yvec];    
                end
                imshow(I);
                hold on;
                plot(pointsx,pointsy,'LineWidth',2);
                plot(x,y,'y.');
                drawnow;
            end
        end
        points=[pointsx; pointsy];
end

% mask calculation
points=points';
if i>2
    mask=roipoly(I,points(:,1),points(:,2));
else
    mask=logical(zeros(siz(1),siz(2)));
end




% internal functions

function [points,t] = myfnplt(f,npoints,varargin)

% interpret the input:
symbol=''; interv=[]; linewidth=[]; jumps=0;
for j=3:nargin
   arg = varargin{j-1};
   if ~isempty(arg)
      if ischar(arg)
         if arg(1)=='j', jumps = 1;
         else, symbol = arg;
         end
      else
         [ignore,d] = size(arg);
         if ignore~=1
	    error('SPLINES:FNPLT:wrongarg',['arg',num2str(j),' is incorrect.']), end
         if d==1
            linewidth = arg;
         else
            interv = arg;
         end
      end
   end
end

% generate the plotting info:
if ~isstruct(f), f = fn2fm(f); end

% convert ND-valued to equivalent vector-valued:
d = fnbrk(f,'dz'); if length(d)>1, f = fnchg(f,'dim',prod(d)); end

switch f.form(1:2)
case 'st'
   if ~isempty(interv), f = stbrk(f,interv);
   else
      interv = stbrk(f,'interv');
   end
%    npoints = 150; 
   d = stbrk(f,'dim');
   switch fnbrk(f,'var')
   case 1
      x = linspace(interv{1}(1),interv{1}(2),npoints);
      v = stval(f,x);
   case 2
      x = {linspace(interv{1}(1),interv{1}(2),npoints), ...
                    linspace(interv{2}(1),interv{2}(2),npoints)};
      [xx,yy] = ndgrid(x{1},x{2});
      v = reshape(stval(f,[xx(:),yy(:)].'),[d,size(xx)]);
   otherwise
      error('SPLINES:FNPLT:atmostbivar', ...
            'Cannot handle st functions with more than 2 variables.')
   end
otherwise
   if ~strcmp(f.form([1 2]),'pp')
      givenform = f.form; f = fn2fm(f,'pp'); basicint = ppbrk(f,'interval');
   end
   
   if ~isempty(interv), f = ppbrk(f,interv); end
      
   [breaks,coefs,l,k,d] = ppbrk(f);
   if iscell(breaks)
      m = length(breaks);
      for i=m:-1:3
         x{i} = (breaks{i}(1)+breaks{i}(end))/2;
      end
%       npoints = 150;
      ii = [1]; if m>1, ii = [2 1]; end
      for i=ii
         x{i}= linspace(breaks{i}(1),breaks{i}(end),npoints);
      end
      v = ppual(f,x);
      if exist('basicint','var') 
                         % we converted from B-form to ppform, hence must now
                         % enforce the basic interval for the underlying spline.
         for i=ii
            temp = find(x{i}<basicint{i}(1)|x{i}>basicint{i}(2));
            if d==1
               if ~isempty(temp), v(:,temp,:) = 0; end
               v = permute(v,[2,1]);
            else
               if ~isempty(temp), v(:,:,temp,:) = 0; end
               v = permute(v,[1,3,2]);
            end
         end
      end
   else     % we are dealing with a univariate spline
%       npoints = 500;
      x = [breaks(2:l) linspace(breaks(1),breaks(l+1),npoints)];
      v = ppual(f,x); 
      if l>1 % make sure of proper treatment at jumps if so required
         if jumps
            tx = breaks(2:l); temp = repmat(NaN, d,l-1);
         else
            tx = []; temp = zeros(d,0);
         end
         x = [breaks(2:l) tx x]; 
         v = [ppual(f,breaks(2:l),'left') temp v];
      end
      [x,inx] = sort(x); v = v(:,inx);
   
      if exist('basicint','var') 
                         % we converted from B-form to ppform, hence must now
                         % enforce the basic interval for the underlying spline.
                         % Note that only the first d components are set to zero
                         % outside the basic interval, i.e., the (d+1)st 
                         % component of a rational spline is left unaltered :-)
         if jumps, extrap = repmat(NaN,d,1); else, extrap = zeros(d,1); end
         temp = find(x<basicint(1)); ltp = length(temp);
         if ltp
            x = [x(temp),basicint([1 1]), x(ltp+1:end)];
            v = [zeros(d,ltp+1),extrap,v(:,ltp+1:end)];
         end
         temp = find(x>basicint(2)); ltp = length(temp);
         if ltp
            x = [x(1:temp(1)-1),basicint([2 2]),x(temp)];
            v = [v(:,1:temp(1)-1),extrap,zeros(d,ltp+1)];
         end
         %   temp = find(x<basicint(1)|x>basicint(2));
         %   if ~isempty(temp), v(temp) = zeros(d,length(temp)); end
      end
   end
   
   if exist('givenform','var')&&givenform(1)=='r' 
                                           % we are dealing with a rational fn:
                                           % need to divide by last component
      d = d-1;
      sizev = size(v); sizev(1) = d;
      % since fnval will replace any zero value of the denominator by 1,
      % so must we here, for consistency:
      v(d+1,find(v(d+1,:)==0)) = 1;
      v = reshape(v(1:d,:)./repmat(v(d+1,:),d,1),sizev);
   end 
end

%  use the plotting info, to plot or else to output:
if nargout==0
   if iscell(x)
      switch d
      case 1
         [yy,xx] = meshgrid(x{2},x{1});
         surf(xx,yy,reshape(v,length(x{1}),length(x{2})))
      case 2
         v = squeeze(v); roughp = 1+(npoints-1)/5;
         vv = reshape(cat(1,...
              permute(v(:,1:5:npoints,:),[3,2,1]),...
              repmat(NaN,[1,roughp,2]),...
              permute(v(:,:,1:5:npoints),[2,3,1]),...
              repmat(NaN,[1,roughp,2])), ...
                    [2*roughp*(npoints+1),2]);
         plot(vv(:,1),vv(:,2))
      case 3
         v = permute(reshape(v,[3,length(x{1}),length(x{2})]),[2 3 1]);
         surf(v(:,:,1),v(:,:,2),v(:,:,3))
      otherwise
      end
   else
      if isempty(symbol), symbol = '-'; end
      if isempty(linewidth), linewidth = 2; end
      switch d
      case 1, plot(x,v,symbol,'linew',linewidth)
      case 2, plot(v(1,:),v(2,:),symbol,'linew',linewidth)
      otherwise
         plot3(v(1,:),v(2,:),v(3,:),symbol,'linew',linewidth)
      end
   end
else
   if iscell(x)
      switch d
      case 1
         [yy,xx] = meshgrid(x{2},x{1});
         points = {xx,yy,reshape(v,length(x{1}),length(x{2}))};
      case 2
         [yy,xx] = meshgrid(x{2},x{1});
         points = {xx,yy,reshape(v,[2,length(x{1}),length(x{2})])};
      case 3
         points = {squeeze(v(1,:)),squeeze(v(2,:)),squeeze(v(3,:))};
         t = {x{1:2}}
      otherwise
      end
   else
      if d==1, points = [x;v];
      else, t = x; points = v([1:min([d,3])],:); end
   end
end


function [xvec,yvec]=EvaluateCardinal2DAtNplusOneValues(P0,P1,P2,P3,T,N)
% Evaluate cardinal spline at N+1 values for given four points and tesion.
% Uniform parameterization is used.

% P0,P1,P2 and P3 are given four points.
% T is tension.
% N is number of intervals (spline is evaluted at N+1 values).

xvec=[]; yvec=[];

% u vareis b/w 0 and 1.
% at u=0 cardinal spline reduces to P1.
% at u=1 cardinal spline reduces to P2.

u=0;
[xvec(1),yvec(1)] =EvaluateCardinal2D(P0,P1,P2,P3,T,u);
du=1/N;
for k=1:N
    u=k*du;
    [xvec(k+1),yvec(k+1)] =EvaluateCardinal2D(P0,P1,P2,P3,T,u);
end


function [xt,yt] =EvaluateCardinal2D(P0,P1,P2,P3,T,u)
% Evaluates 2D Cardinal Spline at parameter value u

% INPUT
% P0,P1,P2,P3  are given four points. Each have x and y values.
% P1 and P2 are endpoints of curve.
% P0 and P3 are used to calculate the slope of the endpoints (i.e slope of P1 and
% P2).
% T is tension (T=0 for Catmull-Rom type)
% u is parameter at which spline is evaluated

% OUTPUT
% cardinal spline evaluated values xt,yt,zt at parameter value u

s= (1-T)./2;

% MC is cardinal matrix
MC=[-s     2-s   s-2        s;
    2.*s   s-3   3-(2.*s)   -s;
    -s     0     s          0;
    0      1     0          0];

GHx=[P0(1);   P1(1);   P2(1);   P3(1)];
GHy=[P0(2);   P1(2);   P2(2);   P3(2)];

U=[u.^3    u.^2    u    1];


xt=U*MC*GHx;
yt=U*MC*GHy;



Contact us