function detbspln
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determines the classification of b-splines for 2D
% n                  : degree
% xx,yy, x,y(gridded): position of grid cells
% D                  : all b-splines supporting the region 
% bb                 : inner(2) and outer(1) b-splines
% gi                 : inner b-splines
% gj                 : outer b-splines
% Created by G.Apaydin in June 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global xx yy n C D bb gi gj x y bound NC;
[x,y]=meshgrid(xx,yy);
if bound{1}.cont==1,e=bound{1}.fun(x,y)>=0;elseif bound{1}.cont==0,e=(x>=min(C(1,:)))&(x<=max(C(1,:)))&(y>=min(C(2,:)))&(y<=max(C(2,:))); else e=isinpoly(x,y,C(1,:),C(2,:));e(e>0)=1;end
if NC>1,for i=2:NC,if bound{i}.cont==1,e=e.*(bound{i}.fun(x,y)>=0);elseif bound{i}.cont==0,e=e.*(x<=min(C(2*i-1,:)))|(x>=max(C(2*i-1,:)))|(y>=min(C(2*i,:)))|(y>=max(C(2*i,:)));else e=e.*(1-isinpoly(x,y,C(2*i-1,:),C(2*i,:)));e(e>0)=1;end,end,end
[mm,nn]=size(x);
%----------------------------DETERMINE BSPLINES----------------------------
bb=zeros([mm,nn]-(n+1));%bbo=bb;bbi=bb;
for i=1:mm-(n+1),for j=1:nn-(n+1),bb(i,j)=sum(sum(e(i:i+n+1,j:j+n+1)))>0;end,end,
% inner bsplines are determined inside all bsplines,
% at the end the remainings are outer bsplines
for i=n+1:mm-(n+1),for j=n+1:nn-(n+1),if (sum(sum(e(i:i+1,j:j+1)))>3.75), bb(i-n:i,j-n:j)=2;end,end,end

[a1,a2]=meshgrid(1:nn-(n+1),1:mm-(n+1));D=[a1(:)';a2(:)';xx(a1(:));yy(a2(:))]';
gi=find(bb==2);gj=find(bb==1);%inner and outer b-splines