Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Update - Thanks everyone for your responses! After fixing a problem with vector normalization, the code below now works.


I'm a new user, and I was attempting to port some Mathematica code from MATLAB for calculating the minimum distance between a point in 3-space and a triangle in 3-space. This is partly a result of me wanting to have this functionality in Mathematica, and, I suppose, a good opportunity to practice learning the Mathematica syntax.

The original MATLAB code isn't mine (and can be found here), but appears to work well in MATLAB. Unfortunately, my translated version of the script in Mathematica returns nonsense despite my almost verbatim copy of the script.

I hope this question isn't inappropriate, but can anyone spot any naïve errors during the translation? Note that I had to change D to Dnp to avoid symbol protection issues.

P = {0.5, -0.3, 0.5};

P1 = {0, -1, 0};
P2 = {1, 0, 0};
P3 = {0, 0, 0};

vertices = {P1, P2, P3};

B = P1;
E0 = P2 - B;
E1 = P3 - B;

Dnp = B - P;
a = Dot[E0, E0];
b = Dot[E0, E1];
c = Dot[E1, E1];
d = Dot[E0, Dnp];
e = Dot[E1, Dnp];
f = Dot[Dnp, Dnp];

det = a*c - b*b;
s = b*e - c*d;
t = b*d - a*e;


If[(s + t) <= det,
  If[s < 0,
    If[t < 0,
      If[(d < 0),
        t = 0;
        If[(-d >= a),
         s = 1;
         sqrDistance = a + 2*d + f;
         ,
         s = -d/a;
         sqrDistance = d*s + f;
         ];
        ,
        s = 0;
        If[(e >= 0),
         t = 0;
         sqrDistance = f;
         ,
         If[(-e >= c),
           t = 1;
           sqrDistance = c + 2*e + f;
           ,
           t = -e/c;
           sqrDistance = e*t + f;
           ];
         ];
        ];
      ,
      s = 0;
      If[e >= 0,
       t = 0;
       sqrDistance = f;
       ,
       If[-e >= c,
         t = 1;
         sqrDistance = c + 2*e + f;
         ,
         t = -e/c;
         sqrDistance = e*t + f;
         ];
       ];
      ];
    ,
    If[t < 0,
      t = 0;
      If[d >= 0,
       s = 0;
       sqrDistance = f;
       ,
       If[-d >= a, 
         s = 1;
         sqrDistance = a + 2*d + f;
         ,
         s = -d/a;
         sqrDistance = d*s + f;
         ];
       ];
      ,
      invDet = 1/det;
      s = s*invDet;
      t = t*invDet;
      sqrDistance = s*(a*s + b*t + 2*d) + t*(b*s + c*t + 2*e) + f;
      ];
    ];
  ,
  If[s < 0,
    tmp0 = b + d;
    tmp1 = c + e;
    If[tmp1 > tmp0,
     numer = tmp1 - tmp0;
     denom = a - 2*b + c;
     If[numer >= denom,
      s = 1;
      t = 0;
      sqrDistance = a + 2*d + f;
      ,
      s = numer/denom;
      t = 1 - s;
      sqrDistance = s*(a*s + b*t + 2*d) + t*(b*s + c*t + 2*e) + f;
      ];
     ,
     s = 0;
     If[tmp1 <= 0,
      t = 1;
      sqrDistance = c + 2*e + f;
      ,
      If[e >= 0,
        t = 0;
        sqrDistance = f;
        ,
        t = -e/c;
        sqrDistance = e*t + f;
        ];
      ];
     ];
    ,
    If[t < 0,
      tmp0 = b + e;
      tmp1 = a + d;
      If[(tmp1 > tmp0),
       numer = tmp1 - tmp0;
       denom = a - 2*b + c;
       If[(numer >= denom),
        t = 1;
        s = 0;
        sqrDistance = c + 2*e + f;
        ,
        t = numer/denom;
        s = 1 - t;
        sqrDistance = s*(a*s + b*t + 2*d) + t*(b*s + c*t + 2*e) + f;
        ];
       ,
       t = 0;
       If[(tmp1 <= 0),
        s = 1;
        sqrDistance = a + 2*d + f;
        ,
        If[(d >= 0),
          s = 0;
          sqrDistance = f;
          ,
          s = -d/a;
          sqrDistance = d*s + f;
          ];
        ];
       ];
      ];
    ,
    numer = c + e - b - d;
    If[numer <= 0,
     s = 0;
     t = 1;
     sqrDistance = c + 2*e + f;
     ,
     denom = a - 2*b + c;
     If[numer >= denom,
      s = 1;
      t = 0;
      sqrDistance = a + 2*d + f;
      ,
      s = numer/denom;
      t = 1 - s;
      sqrDistance = s*(a*s + b*t + 2*d) + t*(b*s + c*t + 2*e) + f;
      ];
     ];
    ];
  ];

If[(sqrDistance < 0), sqrDistance = 0;];

dist = Sqrt[(sqrDistance)]
share|improve this question
    
It's always a good idea to use lower-case letters for your variable names: C, D, E, I, K, N, and O are all reserved. –  cormullion May 16 '13 at 9:54
    
@cormullion Thanks. :) I was just trying to stay "as true" to the original script as possible to make identification of errors easier. So far there haven't been other flags for reserved variable names. –  Richard May 16 '13 at 9:55
1  
It might interest you to know that there is an (undocumented) built-in function, Graphics`Mesh`PointPolygonDistance[] that you might be able to use... –  J. M. May 16 '13 at 11:02
1  
Thanks! - although it's definitely worth waiting a day or two for the real experts to contribute additional answers (I failed maths at school..). As to your original code, I would say (without knowing anything about MatLab) that it's not a good idea to carry out low-level literal translations, but try to find equivalents for the higher-level tasks (vectors, etc). And that nested If statement would look bad in any language... You could usefully read this question too. –  cormullion May 16 '13 at 12:16
2  
@cormullion I don't know about you, but I've never had much success translating procedural code: the abstract problem (if not fundamentally procedural in nature) is often so obscured by the nested loops and conditionals that I have difficulty understanding what is actually going on. There is a tendency in papers for the authors to present algorithms procedurally using buggy, ill-defined pseudocode, which one assumes basically paraphrases their actual program. What is really required in these cases is a flowchart, but these seem to be rather passé nowadays. –  Oleksandr R. May 17 '13 at 3:37

2 Answers 2

up vote 5 down vote accepted

Sorry I can't make sense of your Matlab code - too many nested Ifs. I don't know whether I've adapted belisarius' answer correctly:

point = RandomInteger[100, {3}];
triangle = RandomInteger[100, {3, 3}] ;
lines = Subsets[triangle, {2}];
nline[{start_, end_}, pt_] := 
  Module[{param = ((pt - start).(end - start))/Norm[end - start]^2}, 
    N@{pt, start + Clip[param, {0, 1}] (end - start)}]

Then:

EuclideanDistance @@ nline[#, point] & /@ lines

{67.0659, 78.7704, 67.0103}

Graphics3D[{
 Sphere[point], 
 Polygon[triangle], 
 Red,
 Line[nline[#, point] & /@ lines]}, 
 Axes -> True]

picture

If I've done it right, upvote his answer!, if not, I'll delete mine...

share|improve this answer
    
I'm looking to calculate the minimum distance to the face of the triangle, but still this is very interesting! Thanks! –  Richard May 16 '13 at 11:53

I'm posting this just to show that in Mathematica, it is entirely possible to start from the definitions:

pointTriangleDistance3D[pt_?VectorQ, tri_?MatrixQ] := Module[{p1 = First[tri], s, t}, 
  Sqrt[MinValue[{SquaredEuclideanDistance[pt, p1 + {s, t}.Map[# - p1 &, Rest[tri]]],
                0 <= s <= 1, 0 <= t <= 1, s + t <= 1}, {s, t}]]]

There certainly are more efficient approaches, of course.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.