How to calculate the shortest distance between a point and a line? »

Some problems are easy to solve using pen and paper, but suddenly become hard when you have to express things in mathematics. In this blog we will discuss one such case:

What is the shortest distance between a point and a line (or line segment, or ray)? And what are the coordinates of the point on the line that is closest?

We will go from idea to code in four steps:

1. We sketch the solution for a line, a line segment and a ray.
2. We use linear algebra to define the vector equations that solve the problem in an N-dimensional space.
3. We implement the equations in M-code
4. We turn the code snippets into a robust MATLAB function

Drawing the solution

We solve the problem graphically, by drawing it in 2D. Two points, A and B, define the line, line segment or ray. P is the independent point. The figure below shows how we construct the distance (red, dashed) and closest point (C). We consider three cases:

1. A line through A and B (the top row)
2. A line segment AB (the middle row)
3. A ray which starts at A and passes through B (the bottom row)

The difference between these cases becomes clear when we consider the intersection point of the line through A and B and the perpendicular through P. The closest point can be the intersection point, or A, or B:

1. For a line, the intersection point is the closest point.
2. For a line segment, if the intersection point is on line segment AB, this is the closest point (middle column). If the intersection point is outside segment AB, the closest point is the end of the segment (left and right column).
3. For a ray, if the intersection point is on the ray, this is the closest point (middle and right column). If the intersection point is not on the ray (behind source A), the closest point is source A (left column).

This approach provides insight, but is limited: we need a ruler to measure the distance, an it only works in 2D.

Putting it into math

The next step is to define a set of vector equations that solve the problem in N dimensions.

Let $$\vec{A}$$ and $$\vec{B}$$ be two points on line $$\vec{L}$$. $$\vec{P}$$ is the independent point, which is not on the line.

The vector equation of the line is:

$$\vec{L}(t) = \vec{A} + t\vec{M}$$ (Eq. 1)

with $$t$$ the running parameter and $$\vec{M}$$ the direction vector:

$$\vec{M} = \vec{B} – \vec{A}$$

The running parameter can take different values for a line, line segment and ray:

• A line is infinite,  so $$t \in (-\infty, \infty)$$
• A line segment is finite, so $$t \in \left[0, 1\right]$$
• A ray is semi-infinite, so $$t \in [0, \infty)$$

We need to find the orthogonal projection of our independent point $$\vec{P}$$ on line $$\vec{L}$$. This can be done by using dot products, as follows:

$$\vec{L}(t_0) = \vec{A} + \left[\frac{(\vec{P}-\vec{A}) \cdot \vec{M}}{\vec{M} \cdot \vec{M}}\right] \vec{M}$$ (Eq. 2)

From Eq. 1 and Eq. 2 follows that $$t_0$$, the value of the running parameter at the intersection of $$\vec{L}$$ and its perpendicular through $$\vec{P}$$ is:

$$t_0 = \frac{(\vec{P}-\vec{A}) \cdot \vec{M}}{\vec{M} \cdot \vec{M}}$$

Now that we have all the parameters, we can calculate the shortest distance $$d$$.

In case of a line:

$$d = \left|\vec{P} – (\vec{A} + t_0\vec{M})\right|$$

In case of a line segment:

$$d = \left\{ \begin{array}{ll} \left| \vec{P} – \vec{A} \right| & t_0 \leq 0 \\ \left| \vec{P} – (\vec{A} + t_0\vec{M})\right| & 0 < t_0 < 1 \\ \left| \vec{P} – \vec{B} \right| & t_0 \geq 1 \end{array}\right.$$

And in case of a ray:

$$d = \left\{ \begin{array}{ll} \left| \vec{P} – \vec{A} \right| & t_0 \leq 0 \\ \left| \vec{P} – (\vec{A} + t_0\vec{M})\right| & t_0 > 0 \end{array}\right.$$

Drafting the implementation

We only need a few lines of code to implement the linear algebra equations MATLAB.

• Direction vector $$\vec{M}$$:
• Running parameter at the orthogonal intersection $$t_0$$:
• Intersection point:
• Distance $$d$$ between independent point and closest point:

Implementing a function

The last step involves coding a robust, documented, and readable MATLAB function. At least three input arguments are required: the points A and B that define a line and test point P. The optional fourth input argument specifies the line type: 'line' (the default), 'segment' or 'ray'. The function returns up to three outputs: distance d, closest point C, and running parameter at the orthogonal intersection t0. The function definition line reflects this:

Download the distancePoint2Line() function here to inspect the full implementation and to use it for your work.

Simulink and Stateflow offer ways to customize the modeling environment. You may for example add a menu to the Simulink Model Editor, or add menu items to existing Model Editor menus. Customizations are a practical way to create shortcuts, allowing you to be more productive in Simulink. However, it can be hard to keep a clear oversight of the customization implementations. The two main reasons are:

1. Each customization needs a function to implement the callback, and a function to register it. In addition, separate functions are required to define and register the menu with the Simulink Customization Manager. Consequently, if you have multiple customizations for diffent projects or contexts, the number of functions will grow rapidly, and significant bookkeeping is required.
2. If you work on different projects you might want to make a selection of the customizations you find useful. Or, you might want to target customizations to a specific context. In these cases, the default approach, a single sl_customization.m file with numerous sub functions, will be cumbersome. The alternative approach, multiple sl_customization.m files in different folders, may lead to a messy folder structure.

This calls for a structured approach. Below we present one, and share the code.

EDIT 22 March 2017: Code and MATLAB App can be found on GitHub and MATLAB Central File Exchange from now.

Example

The first, “Match Size”, matches the size of selected blocks. This can also be arranged via the Simulink context menu, however the customization adds a keyboard shortcut to it.

The second, “Make Constants Orange”, colors all Constant blocks in a model orange, which reflects a convention.

The unstructured approach

All code is in one sl_customization.m file, which has six sub functions. It works, but if you want to add more customizations this file will soon grow to a size that is no longer maintainable.

The structured approach

The two customizations in our example originate from different contexts: one is a general utility, the other part of a style guide. In the structured approach, separate packages (+styleguide, +utils) reflect the different contexts. Both packages contain a Customizer class. A single sl_customization.m file points to the Customizers employed in our model.

The idea is that we separate the definition of the custom menu and its contents from the implementation of the customizations. This way, we can add or remove menu items by editing a single line in sl_customization.m. Also, if we add customizations to a package, they will automatically appear in our custom menu.

The implementation reflects this separation. A generic class Customizer is the super class of the Customizer classes in the package folders. The super class contains a helper method getCustomizeMethods that will automatically assign the customizations to the menu. The customizers variable in sub-function customMenu (line 27) defines the items in menu “My Menu”. The Customizer super class will place a separator between customizations from different packages.

The Customizer class in the +utils package implements the details for the “Match Size” menu item: a registration method matchSizeOfBlocks and a callback sub-function matchSizeOfBlocksCb.

Similarly, the Customizer class in the +styleguide package implements the details for the “Make Constants Orange” menu item.

In short, this structure allows to easily scale up the number of customizations. Also, it facilitates sharing and re-use. The full implementation of the structured approach is available for download. Use it to automate your repeat work and be more productive in Simulink/Stateflow!

EDIT 22 March 2017: Code and MATLAB App can be found on GitHub and MATLAB Central File Exchange from now.

MonkeyProof Solutions @ MATLAB Expo 2016 Benelux »

MonkeyProof Solutions exhibits at MATLAB Expo 2016 Benelux in Eindhoven on Tuesday June 28, 2016.

Meet us at our booth! This is an excellent opportunity to discuss your MATLAB, modeling, and Model-Based Design challenges and/or needs with the experts of MonkeyProof Solutions.

We invite you to send us a quick message ahead of the event such that we can make sure that you meet the right expert!

Speeding up a MATLAB simulation – A case study »

The benefit of speed

In model-based design, the model is the workhorse for product development. In general, the main objective during model definition is to implement the required features and functions. Speed of execution may be a lower priority. In practice, improving run-time performance is often postponed, as it takes time and may require re-programming the fresh implementation. However, it may be worthwhile: short execution times are beneficial both in the development and the deployment phase.

During development, a short execution time provides timely feedback: it invites people to test more implementation ideas, and to more thoroughly validate the model implementation. The resulting quick iterations will lead to higher quality models.

During deployment, one can gain much more insight from a speedy model: explore numerous what-if scenarios, explore parameter spaces, optimize parameters, run more complex cases, validate people’s gut feeling during discussions, etc. A high-performing model may even become part of real-time applications, or run on a server to be deployed by many.

In this blog we provide a real-world example of model speedup: we outline a possible approach and  show how simulation time was accelerated 15 times, while making no concessions to simulation accuracy.

Case: A subsidence model

In a recent project, we implemented a subsidence model. Subsidence is the vertical movement of the ground surface due to pressure reduction in a reservoir. Subsidence can be a side-effect of gas exploitation. Van Opstal’s method is an often used, semi-analytical approach to predict subsidence.  Van Opstal’s equations quantify the contribution of a reservoir cell at position (x,y,z) to the eventual subsidence of a point (x,y) at the surface. In our model, the cells that compose an arbitrary 3D reservoir shape have an in-homogeneous  and time-dependent pressure distribution. At the ground surface we define a rectangular grid around the reservoir footprint. We use the model to study subsidence over time, by feeding it time-varying reservoir cell pressures P(x,y,z,t).

M-code implementation of Van Opstal’s equations

We implemented the subsidence model in a few MATLAB classes and translated Van Opstal’s equation set to a solver method. Our base case contains 20k surface grid elements and 115k reservoir elements. As a result, a total 2.3G subsidence contributions need to be calculated to run the case. On the laptop (Intel Core2 Duo T9400 @2.53GHz) the simulation took 163 minutes (2h43m) to complete.

Speedup options

There are many options to speed up MATLAB code in general. Which options are effective depends on factors like functions used, data dimensions, numerical accuracy, MATLAB version, hardware, bottlenecks related to I/O, etc. It is even possible to run code on a computer cluster using the MATLAB Distributed Computing Server.

In this blog we focus on a single box solution, and apply and quantify four options to speed up the model:

1. Faster hardware
2. M-code optimization
3. Parallel processing on CPU
4. Parallel processing on GPU

Option 1: Faster hardware

A logical first step is to run the simulation on faster hardware. We ran the model on a quad-core machine (Intel Core i7 CPU @2.93 GHz, 4 physical cores — an average CPU). The simulation now took 45 minutes, a 3.6x speedup. This will be our performance benchmark for further work.

Option 2: M-code optimization

We used the MATLAB profiler to identify performance bottlenecks. The main bottleneck appeared to be the calculation of factors that include element-wise divisions and element-wise powers. Especially the element-wise power is expensive in MATLAB, because MathWorks uses the fdlibm library which is more accurate but slower than straight multiplication (see also Bobby Cheng’s post). We re-coded some sections and validated the accuracy of the implementation. The simulation now took 8 minutes, a further 5.6x speedup.

Option 3: Parallel processing on CPU

The Van Opstal solver spends most time in a for loop which calculates subsidence at each surface grid point. The iterations of this loop are independent which means that the problem is massively parallel. We can turn the for loop in a parfor loop, and have a so-called parallel pool of workers (multiple headless MATLAB sessions running on separate cores) perform the iterations, using all cores of a machine.

However, it should be taken into account that recent MATLAB versions implement implicit multi threading for many base functions. In our case, vectors of more than 100k elements (that is: the reservoir data vectors) trigger multi threaded behavior for operations like addition, subtraction, multiplication, division and square root.  Therefore, the parfor construct will not speed up our code on a single machine. Actually, a parallel pool of workers slowed down the computations due to the overhead it introduces. We better skip this option.

Option 4: Parallel processing on GPU

General purpose Graphics Processing Units (GPUs) are high-performance many-core processors that can be used to accelerate engineering computations. We ran our simulation on a powerful NVIDIA graphics cards using the GPU interface provided by the Parallel Computing Toolbox.

We systematically assessed the GPU acceleration potential of individual operations and functions to estimate the overall GPU acceleration potential for the Van Opstal solver. In this case we expected an overall speedup of 10-15x for a mid-range NVIDIA Quadro GPU card.

However, it is not always straightforward to release the full GPU potential from MATLAB. Simply transferring the MATLAB arrays to the GPU (using the gpuArray class) and letting it do the math does not work in this case. On the contrary: that approach slowed down the solver.

Instead, we modified our solver implementation in order to feed the GPU tailored chunks. We turned the vectorized implementation into an element-wise implementation, and used gpuArray.arrayfun() to calculate a number of array elements in parallel on the GPU. After tuning the number of parallel GPU threads, it ran stable and fast, resulting in a 14.9x acceleration compared to benchmark. One model run takes 3 minutes to complete. This is acceptable performance: we can now simulate a 30-year time horizon in 1.5 hours instead of a day.

Summary

In short, a GPU implementation accelerated our model 14.9 times compared to the original M-code running on a quad-core machine, and 53.7 times compared to running it on a dual-core laptop. It was also shown that there is no “one size fits all” approach: potential speedup options may prove ineffective, or may require tuning.

Speed increases the quality, applicability and therefore value of a model. Do you need to speed up or enhance your simulations? Contact MonkeyProof Solutions to discuss how the experts can be of help.