Index: /issm/trunk/src/m/classes/public/processgeometry.m
===================================================================
--- /issm/trunk/src/m/classes/public/processgeometry.m	(revision 3665)
+++ /issm/trunk/src/m/classes/public/processgeometry.m	(revision 3666)
@@ -14,6 +14,7 @@
 	x2=geom.Vertices(geom.Edges(i,2),1);
 	y2=geom.Vertices(geom.Edges(i,2),2);
+	color1=geom.Edges(i,3);
 
-	j=i+1; %test edges located AFTER i only
+	j=i; %test edges located AFTER i only
 	while (j<size(geom.Edges,1)),
 
@@ -31,4 +32,5 @@
 		x4=geom.Vertices(geom.Edges(j,2),1);
 		y4=geom.Vertices(geom.Edges(j,2),2);
+		color2=geom.Edges(j,3);
 
 		%Check if the two edges are crossing one another
@@ -40,5 +42,5 @@
 
 			%Add vertex to the list of vertices
-			geom.Vertices(end+1,:)=[x y 1];
+			geom.Vertices(end+1,:)=[x y min(color1,color2)];
 			id=size(geom.Vertices,1);
 
@@ -56,4 +58,39 @@
 	end
 
+end
+
+%Check point outside
+disp('Checking for points outside the domain...');
+i=0;
+num=0;
+while (i<size(geom.Vertices,1)),
+
+	%vertex counter
+	i=i+1;
+
+	%Get coordinates
+	x=geom.Vertices(i,1);
+	y=geom.Vertices(i,2);
+	color=geom.Vertices(i,3);
+
+	%Check that the point is inside the domain
+	if (color~=1 & ~ContourToNodes(x,y,outline(1),1)),
+
+		%Remove points from list of Vertices
+		num=num+1;
+		geom.Vertices(i,:)=[];
+
+		%update edges
+		[posedges dummy]=find(geom.Edges==i);
+		geom.Edges(posedges,:)=[];
+		posedges=find(geom.Edges>i);
+		geom.Edges(posedges)=geom.Edges(posedges)-1;
+
+		%update counter
+		i=i-1;
+	end
+end
+if num,
+	disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']);
 end
 
@@ -102,36 +139,2 @@
 %remove empty edges
 geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
-
-%Check point outside
-disp('Checking for points outside the domain...');
-i=0;
-num=0;
-while (i<size(geom.Vertices,1)),
-
-	%vertex counter
-	i=i+1;
-
-	%Get coordinates
-	x=geom.Vertices(i,1);
-	y=geom.Vertices(i,2);
-
-	%Check that the point is inside the domain
-	if (~ContourToNodes(x,y,outline(1),1)),
-
-		%Remove points from list of Vertices
-		num=num+1;
-		geom.Vertices(i,:)=[];
-
-		%update edges
-		[posedges dummy]=find(geom.Edges==i);
-		geom.Edges(posedges,:)=[];
-		posedges=find(geom.Edges>i);
-		geom.Edges(posedges)=geom.Edges(posedges)-1;
-
-		%update counter
-		i=i-1;
-	end
-end
-if num,
-	disp(['WARNING: ' num2str(num) 'points outside the domain outline have been removed']);
-end