13
$\begingroup$

I'm using version 11 for whatever that's worth. I'm doing an optical simulation with thousands of converging rays (HalfLines), since this is in 3D for the most part they don't intersect, but I would like to find the closest point among all of the rays. I'm new to Mathematica so I'm trying to figure out what an efficient and idiomatic way to do this is. The docs have a nice example with two regions like this:

{d, pqrul} = Minimize[EuclideanDistance[p, q], {p ∈ 𝒫, q ∈ 𝒬 // RootReduce

I've thought of trying to adapt that somehow (all pairs of points in the RegionUnion of the lines?) but that seems grossly inefficient. There's always the MATLAB way. Still, I'd like to learn the idiomatic way because I'd also like to be able to visualize these things without packing them into and stripping them out of matrices.

Improve this question
$\endgroup$
4
  • $\begingroup$ Have you seen RegionDistance[] and RegionNearest[]? $\endgroup$
    – J. M.'s eventual burnout ♦
    Aug 13, 2016 at 2:37
  • 1
    $\begingroup$ I already wrote an answer, but became unsure afterwards: do you want your solution to be on some of the HalfLines, or just the least squares solution to the problem, anywhere in space? $\endgroup$
    – kirma
    Aug 13, 2016 at 6:02
  • $\begingroup$ It can be anywhere in space, your answer looks great! Thanks! $\endgroup$
    – Dan
    Aug 13, 2016 at 6:27
  • 1
    $\begingroup$ Do you need HalfLines, or would InfiniteLines work, too? That would simplify the problem a lot. I ask because the "Matlab way", as far as I can see, assumes infinite lines. $\endgroup$
    – Niki Estner
    Aug 13, 2016 at 7:05

2 Answers 2

Reset to default
13
$\begingroup$

EDIT: As @nikie noted, using FindArgMin (a variant of FindMinimum) instead of (N)ArgMin can improve the speed of finding a solution. Since in the case of this problem only one minimum exists, this should actually produce the global minimum in an efficient manner.

How about this?

Module[{n, lines, sol},
 n = 10;
 lines = HalfLine[# + RandomPoint@Ball[], -2 # + RandomPoint@Ball[]] & /@
      RandomPoint[Sphere[{0, 0, 0}, 10], n];
 sol = 
  FindArgMin[Norm[RegionDistance[#, {x, y, z}] & /@ lines], {x, y, z}];
 Graphics3D[{Arrow@Tube@{#1, #1 + #2} & @@@ lines, 
   Red, PointSize@Large, Point@sol}]]

enter image description here

Here lines is a list of random almost-focused HalfLines, sol is the coordinate of least-squares solution, and rest is just visualizing the result (red Point marks the spot).

This solution is constrained - by the fact HalfLines have a starting point, and that distances behind this point are measured to it, instead of the projection on the infinite line. This differs from solution provided by @J.M. which employs fully infinite lines. My solution can be turned into equivalent with that solution just by replacing HalfLine with InfiniteLine.

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ Er, I'm not sure why, but it seems that replacing ArgMin into {x,y,z}/.NMinimize[***][[2]] can increase the speed for 3 times when 30 lines are presented. $\endgroup$
    – Wjx
    Aug 13, 2016 at 6:40
  • 1
    $\begingroup$ Oh, yeah yeah, I didn't see that there's no N before ArgMin. I suppose this sort of things won't happen? The documentation said that NArgMin will always try to find the global minimum. $\endgroup$
    – Wjx
    Aug 13, 2016 at 6:47
  • 1
    $\begingroup$ FindMinimum seems much faster than NArgMin. It displays a warning message, but returns the same result. And if you could replace HalfLines with InfiniteLines, it gets much faster, but I'm not sure if that's required for the question. $\endgroup$
    – Niki Estner
    Aug 13, 2016 at 7:07
  • 3
    $\begingroup$ @kirma: But there is only one minimum, right? FindMinimum is guaranteed to find a minimum, so if there is only one - how could FindMinimum not find it? $\endgroup$
    – Niki Estner
    Aug 13, 2016 at 7:20
  • 2
    $\begingroup$ @MichaelE2: You're right, I was speaking imprecisely: The objective function here is convex (the euclidean distance from a [half-]line is convex, and the sum of convex functions is also convex). FindMinimum should (in theory) always find the single, global minimum for this class of objective. In practice, if you're including terms like x^100, it might not, due to numerical overflow or underflow. But you certainly could come up with equivalent examples for NMinimize $\endgroup$
    – Niki Estner
    Aug 13, 2016 at 13:48
11
$\begingroup$

Here's another way to implement kirma's solution:

With[{n = 15}, BlockRandom[SeedRandom["many lines"]; (* for reproducibility *)
     lines = HalfLine[# + RandomPoint @ Ball[], -2 # + RandomPoint @ Ball[]] & /@ 
     RandomPoint[Sphere[{0, 0, 0}, 10], n]]];

sol = LeastSquares @@ Total[
      With[{m = IdentityMatrix[3] - (Outer[Times, #, #] &[Normalize[#2]])},
           {m, m.#1}] & @@@ lines];

Graphics3D[{{Gray, Arrow[Tube[{#1, #1 + #2} & @@@ lines]]},
            {Directive[Red, PointSize[Large]], Point[sol]}}, Boxed -> False]

least-squares "intersection" point

Improve this answer
$\endgroup$
4
  • $\begingroup$ It's not really too obvious if this implementation handles HalfLines as what they are defined to be, or as rays extending on both directions to the infinity. Which way it is? $\endgroup$
    – kirma
    Aug 13, 2016 at 9:09
  • 1
    $\begingroup$ @kirma, that would be the second one; as I said, I just reimplemented your solution in terms of LeastSquares[]. You can check that this returns the same solution as yours. $\endgroup$
    – J. M.'s eventual burnout ♦
    Aug 13, 2016 at 9:22
  • 3
    $\begingroup$ I posted an example of differing results (half-lines vs. "complete" ones) on chat ( i.stack.imgur.com/cx7BE.png ). Our answers are solutions to different questions. :) $\endgroup$
    – kirma
    Aug 13, 2016 at 9:33
  • 1
    $\begingroup$ Thank you for this solution also. I like how nicely it maps on to the Wikipedia notation (en.wikipedia.org/wiki/… ). $\endgroup$
    – Dan
    Aug 13, 2016 at 21:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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