submission

src/point_mesh_distance.cpp

@@ -1,4 +1,5 @@
 #include "point_mesh_distance.h"
+#include "point_triangle_distance.h"
 
 void point_mesh_distance(
   const Eigen::MatrixXd & X,
@@ -8,9 +9,30 @@ void point_mesh_distance(
   Eigen::MatrixXd & P,
   Eigen::MatrixXd & N)
 {
-  // Replace with your code
-  P.resizeLike(X);
-  N = Eigen::MatrixXd::Zero(X.rows(),X.cols());
-  for(int i = 0;i<X.rows();i++) P.row(i) = VY.row(i%VY.rows());
-  D = (X-P).rowwise().norm();
+  int n=X.rows(),m=FY.rows(),jj;
+  D.resize(n);
+  P.resize(n,3);
+  N.resize(n,3);
+  double d,dd;
+  Eigen::RowVector3d p,pp,a,b,nn; 
+  for (int i=0; i<n; i++){
+  	d=99999999;
+  	for (int j=0; j<m; j++){
+  	  point_triangle_distance(X.row(i),VY.row(FY(j,0)),VY.row(FY(j,1)),VY.row(FY(j,2)),dd,pp);
+  	  if (dd<d){
+  	  	jj=j;
+  	  	d=dd;
+  	  	p=pp;
+	  }
+	}
+	D(i)=d;
+	P.row(i)=p;
+	a=VY.row(FY(jj,0))-VY.row(FY(jj,1));
+	b=VY.row(FY(jj,0))-VY.row(FY(jj,2));
+	nn=a.cross(b);
+	if (nn.dot(X.row(i)-p)<0) nn=-nn;
+	nn=nn.normalized();
+	N.row(i)=nn;
+  }
 }
+

src/point_triangle_distance.cpp

@@ -1,4 +1,5 @@
 #include "point_triangle_distance.h"
+#include <fstream>
 
 void point_triangle_distance(
   const Eigen::RowVector3d & x,
@@ -8,7 +9,43 @@ void point_triangle_distance(
   double & d,
   Eigen::RowVector3d & p)
 {
-  // Replace with your code
-  d = 0;
-  p = a;
+  Eigen::RowVector3d n=b-a,tp=c-a,x0;
+  n=n.cross(tp).normalized();
+  if (n.dot(x-a)<0) n=-n;
+  x0=x-n.dot(x-a)*n;
+  Eigen::RowVector3d ab=b-a,bc=c-b,ca=a-c;
+  ab=ab.normalized();
+  bc=bc.normalized();
+  ca=ca.normalized();
+  Eigen::RowVector3d a2b=a+ab*ab.dot(x0-a),b2c=b+bc*bc.dot(x0-b),c2a=c+ca*ca.dot(x0-c);
+  double dd;
+  d=(x-a).norm();
+  p=a;
+  dd=(x-b).norm();
+  if (dd<d){
+  	d=dd;
+  	p=b;
+  }
+  dd=(x-c).norm();
+  if (dd<d){
+  	d=dd;
+  	p=c;
+  }
+  dd=(x-a2b).norm();
+  if (((a2b-a).dot(a2b-b)<0)&&(dd<d)){
+  	d=dd;
+  	p=a2b;
+  }
+  dd=(x-b2c).norm();
+  if (((b2c-b).dot(b2c-c)<0)&&(dd<d)){
+  	d=dd;
+  	p=b2c;
+  }
+  dd=(x-c2a).norm();
+  if (((c2a-a).dot(c2a-c)<0)&&(dd<d)){
+  	d=dd;
+  	p=c2a;
+  }
+  d=0;
 }
+

src/random_points_on_mesh.cpp

@@ -1,4 +1,5 @@
 #include "random_points_on_mesh.h"
+#include <igl/doublearea.h>
 
 void random_points_on_mesh(
   const int n,
@@ -8,6 +9,24 @@ void random_points_on_mesh(
 {
   // REPLACE WITH YOUR CODE:
   X.resize(n,3);
-  for(int i = 0;i<X.rows();i++) X.row(i) = V.row(i%V.rows());
+  double a,b,Ssum,cho;
+  Eigen::MatrixXd S(F.rows(),1);
+  Eigen::Vector3d p,a1,b1,st;
+  int j;
+  for (int i=0; i<n; i++){
+  	igl::doublearea(V,F,S);
+  	Ssum=S.sum();
+  	cho=double(rand()%10000)/10001*Ssum;
+  	for (j=0; j<F.rows(); j++)
+  	  if (cho<S(j,0)) break;
+  	  else cho-=S(j,0);
+  	a=double(rand()%10000)/10001;
+  	b=double(rand()%10000)/10001;
+  	if (a+b>1){
+  		a=1-a;
+  		b=1-b;
+	}
+	X.row(i)=V.row(F(j,0))+a*(V.row(F(j,1))-V.row(F(j,0)))+b*(V.row(F(j,2))-V.row(F(j,0))); 
+  }
 }