function  levelset=exp_to_levelset(md,contourname)
%EXP_TO_LEVELSET - from an exp contour file, build a levelset which has positive sign inside the contour, and negative outside.
%
%   Usage:
%      levelset=exp_to_levelset(md,contourname)
%
%   Where: 
%      - 'md'     : model specifying the mesh
%      - 'contourname'  : exp file which defines the closed contour
%      - 'levelset'  : levelset vector on each vertex of the mesh, size md.mesh.numberofvertices
%
%   Example:
%      level=exp_to_levelset(md,'DomainOutline.exp');
%
%   See also MESHPROFILEINTERSECTION, SHPREAD, EXPREAD
	
	tolerance=1e-5; %tolerance level used for detection of edges, points coinciding, etc ...
		
	%from exp contour, build list of segments throughout our mesh:
	segments=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,contourname);

	%cleanup: remove 0 length segments: 
	cleanup=1;
	if cleanup,
		flags=zeros(size(segments,1),1);
		for j=1:size(segments,1),
			segment=segments(j,:);
			x1=segment(1); x2=segment(3); 
			y1=segment(2); y2=segment(4); 
			Q1=[x1;y1];Q2=[x2;y2];
			if norm(Q2-Q1)==0, 
				flags(j)=1;
			end
		end
		pos=find(flags);
		segments(pos,:)=[];
	end

	%simplify segments. levelsets are P1 through our triangular elements, which means we can only have one segment per element. 
	%We simplify by merging all segments into one for each element.
	simplify=1;
	if simplify,
		ids=unique(segments(:,5));
		newsegs=zeros(length(ids),5);
		newsegs(:,5)=ids;
		for j=1:length(ids),
			pos=find(segments(:,5)==ids(j));
			segs=segments(pos,1:4);
			count=0;
			while(size(segs,1)>1),
				count=count+1;
				if count>10,
					%keep first segment cause something funky in the configuration is happening.
					segs=segs(1,:);
					break;
				end
				if (abs(segs(1,3)-segs(2,1))<tolerance & abs(segs(1,4)-segs(2,2))<tolerance),
					%the second segment disappears, extending the first segment:
					if size(segs,1)==2,
						segs=[segs(1,1) segs(1,2) segs(2,3) segs(2,4)];
					else
						segs=[segs(1,1), segs(1,2), segs(2,3), segs(2,4); segs(3:end,:)];
					end
				else
					%put the first segment at the end
					segs=[segs(2:end,:); segs(1,:)];
				end
			end
			newsegs(j,1:4)=segs;
		end
		segments=newsegs;
	end

	%cleanup: remove 0 length segments: 
	if cleanup,
		flags=zeros(size(segments,1),1);
		for j=1:size(segments,1),
			segment=segments(j,:);
			x1=segment(1); x2=segment(3); 
			y1=segment(2); y2=segment(4); 
			Q1=[x1;y1];Q2=[x2;y2];
			if norm(Q2-Q1)<tolerance,
				flags(j)=1;
			end
		end
		pos=find(flags);
		segments(pos,:)=[];
	end

	%now, using the segments, create level set flags on each element that touches the segment.
	distance=zeros(md.mesh.numberofvertices,1);
	flags=zeros(md.mesh.numberofvertices,1);
	for j=1:size(segments,1),
		segment=segments(j,:);
		id=segment(5);
		x1=segment(1); x2=segment(3); 
		y1=segment(2); y2=segment(4); 
		for k=1:3,
			node=md.mesh.elements(id,k);
			xk=md.mesh.x(node); yk=md.mesh.y(node);
			
			%on which side of the segment is this node k? 
			P=[xk;yk]; Q1=[x1;y1];Q2=[x2;y2];
			d=(det([Q2-Q1,P-Q1]))/norm(Q2-Q1);
			if(~flags(node)),
				distance(node)=d;
				flags(node)=1;
			end
		end
	end
		
	%now go through all the other nodes and complete (similar to phase unwrapping: we march through the domain, and each time we 
	%are connected to a node which has a levelset value, we adopt the same one.
	complete=1;
	if complete,
		while(length(find(flags==0))),
			sumf=sum(flags(md.mesh.elements),2);
			pos=find(sumf==2);
			if ~isempty(pos),
				for j=1:length(pos),
					el=pos(j);
					f=flags(md.mesh.elements(el,:));
					ids=find(f==1); 
					d=distance(md.mesh.elements(el,ids));
					if d(1)*d(2)<0,
						%we have hit a snag here, can't decide!
						continue;
					else
						d=d(1);
					end
					if abs(d)<tolerance,
						continue;
					end
					ids=find(f==0);
					if ~isempty(ids),
						flags(md.mesh.elements(el,ids))=1;
						distance(md.mesh.elements(el,ids))=d;
					end
				end
			end
		end
	end
	levelset=distance;
