Index: /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp	(revision 22455)
+++ /issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp	(revision 22456)
@@ -9,9 +9,8 @@
 #endif
 
-/*Include files: {{{*/
+/*Include files*/
 #include "./ExpToLevelSetx.h"
 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0);
-void ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
-/*}}}*/
+void   ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
 
 void* ExpToLevelSetxt(void* vpthread_handle){
@@ -29,9 +28,9 @@
 
 	/*recover parameters :*/
-	Contours* contours  = gate->contours;
-	int       nods      = gate->nods;
-	double   *distance    = gate->distance;
-	double   *x         = gate->x;
-	double   *y         = gate->y;
+	Contours* contours = gate->contours;
+	int       nods     = gate->nods;
+	double   *distance = gate->distance;
+	double   *x        = gate->x;
+	double   *y        = gate->y;
 
 	/*distribute indices across threads :*/
@@ -39,5 +38,5 @@
 
 	/*Loop through all contours: */
-	for (i=0;i<contours->Size();i++){
+	for(i=0;i<contours->Size();i++){
 		Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
 		ContourToLevelSet(distance,contour->x,contour->y,contour->nods,x,y,i0,i1);
@@ -48,5 +47,4 @@
 
 void ContourToLevelSet(double* dist,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i1){/*{{{*/
-	int i,j;
 	double x0,y0;
 	double x1,y1;
@@ -54,11 +52,14 @@
 	double mind;
 	
-	for(i=i0;i<i1;i++){
+	for(int i=i0;i<i1;i++){
+
+      /*Get current point*/
 		x0=x[i]; y0=y[i];
 
 		/*Figure out distance from (x0,y0) to contour: */
 		mind=1e+50;
-		for (j=0;j<contournods-1;j++){
-			x1=contourx[j]; y1=contoury[j];
+		for(int j=0;j<contournods-1;j++){
+         /*Get distance from current segment*/
+			x1=contourx[j];   y1=contoury[j];
 			x2=contourx[j+1]; y2=contoury[j+1];
 			mind=min(mind,minimum_distance(x1,y1,x2,y2,x0,y0));
@@ -66,28 +67,32 @@
 		dist[i]=min(dist[i],mind);
 	}
-
 }
 double ddistance(double x1,double y1,double x2,double y2){ return sqrt(pow(x2-x1,2)+pow(y2-y1,2)); }
 double ddot(double x1, double y1, double x2, double y2){ return x1*x2+y1*y2; }
 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
+	/* Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)*/
+
+   /*Get segment length square (avoid sqrt) |w-v|^2*/
+	double l2 = pow(x2-x1,2)+pow(y2-y1,2); 
+
+   /*segment is single point: v == w*/
+	if(l2 == 0.) return ddistance(x0,y0,x1,y1);
+
+	/*Consider the line extending the segment, parameterized as v + t (w - v).
+	We find projection of point p onto the line. 
+	It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+	double t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2;
+	if(t < 0.){
+      /*Beyond the 'v' end of the segment*/
+      return ddistance(x0,y0, x1, y1);
+   }
+	else if(t > 1.){
+      /*Beyond the 'w' end of the segment*/
+      return ddistance(x0,y0, x2,y2);
+   }
 	
-	// Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)
-	double projectionx; 
-	double projectiony; 
-	double l2;
-	double t;
-
-	l2 = pow(x2-x1,2)+pow(y2-y1,2); // i.e. |w-v|^2 -  avoid a sqrt
-
-	if (l2 == 0.0) return ddistance(x0,y0, x1,y1); // v == w case
-	// Consider the line extending the segment, parameterized as v + t (w - v).
-	//         // We find projection of point p onto the line. 
-	//           // It falls where t = [(p-v) . (w-v)] / |w-v|^2
-	t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2;
-	if (t < 0.0) return ddistance(x0,y0, x1, y1);       // Beyond the 'v' end of the segment
-	else if (t > 1.0) return ddistance(x0,y0, x2,y2);  // Beyond the 'w' end of the segment
-	
-	projectionx= x1 + t* (x2-x1);  // Projection falls on the segment
-	projectiony= y1 + t* (y2-y1);
+   /*Projection falls on segment!*/
+	double projectionx= x1 + t* (x2-x1);
+	double projectiony= y1 + t* (y2-y1);
 	return ddistance(x0, y0, projectionx, projectiony);
 }