One of the techniques I developed for my thesis research was called parabolic arc-length projection. The purpose of this projection was to help implement a hierarchical image segmentation technique for dual-polarization RADARSAT images. In the course of coding the projection, I learned how to vectorize analytical expressions found by MATLAB's Symbolic Solve. This lets you use the symbolic expressions with vector data, allowing you to find the value of each expression for multiple input values without using loops. Avoiding loops is essential because they tend to be slow in MATLAB.

Given a series of 2D points and the equation of a parabola (`y = a1*x^2 + b1*x + c`

), find for each point the nearest neighbouring point on the parabola. The figure below shows this for one point, `p = (x1, y1)`

:

The nearest point on the parabola is found analytically by minimizing the Euclidean distance equation between point `p`

and a generic point on the parabola whose coordinates are `(x, a1*x^2 + b1*x + c1)`

. If you work out the math, you end up with a general cubic equation to solve, whose coefficients `a`

, `b`

, `c`

and `d`

, are functions of `x1`

, `y1`

, `a1`

, `b1`

and `c1`

. Solving the equation gives `x'`

, the x coordinate of the point on the parabola that is closest to point `p`

. The coordinates of this point are `(x', a1*x'^2 + b1*x' + c1)`

The formulae for solving cubic equations can be found in G. Cardano's *Artis magnae, sive de regulis algebraicis liber unus* (Book one about the great art / the rules of algebra). I spent a bit of time trying to translate *Artis magnae* but it was a bit more challenging to parse than *Cambridge Latin Course*, *De Vita Caesarum* or *Commentarii de Bello Gallico*. Since my Latin was a bit rusty, I ended up using MATLAB's symbolic solver. The Symbolic Solver gives you a symbolic expression that is the solution of the equation you want to solve.

A generic cubic equation in MATLAB can be expressed as a symbolic expression and then solved using `solve`

(you will need the Symbolic Math Toolbox):

% declare symbolic variables. % a, b, c and d are coefficients of the cubic to be solved. syms a b c d x % use the symbolic solver to solve a generic cubic equation. % a*x^3+b*x^2+c*x+d = 0 r=solve(a*x^3+b*x^2+c*x+d);

The variable `r`

is actually an array of three symbolic expressions, corresponding to equations of the three roots of the cubic equation. Take a look:

>> r(1) ans = 1/6/a*(36*c*b*a-108*d*a^2-8*b^3+12*3^(1/2)*(4*c^3*a-c^2*b^2-18*c*b*a*d+27*d^2*a^2+4*d*b^3)^(1/2)*a)^(1/3)- [...]

I truncated the expression because it's very long.

To find the roots of the cubic equation, all you need to do is set the values of `a`

, `b`

, `c`

and `d`

(which are functions of the point you are trying to analyze and the equation coefficients of the parabola) and then compute the value of each symbolic expression with this set of `a`

, `b`

, `c`

and `d`

values.

To use the symbolic expression as a function, create inline functions using the symbolic expressions by using `inline`

:

R1_scalar = inline(char(r(1)),'a','b','c','d'); R2_scalar = inline(char(r(2)),'a','b','c','d'); R3_scalar = inline(char(r(3)),'a','b','c','d');

Each expression (e.g. `r(1)`

) is converted to a string first using `char`

. This puts it in the right format for `inline`

to create a function out of it. The `inline`

calls create functions `RN_scalar`

that takes four arguments, `(a, b, c, d)`

. Once you have set the values for `a`

, `b`

, `c`

, `d`

for the 2D point for which you are trying to find the nearest parabola point, calling `RN_scalar(a,b,c,d)`

will give you the corresponding root.

For the problem at hand - implementing the parabolic arc-length projection algorithm - there are many points for which the nearest neighbour on a parabolic curve must be found. Since this is an image processing application, there is one 2D point per pixel of an image, which means that there are *many* points that need to have the equation applied:

Point | Coordinates | Cubic Coefficients |
---|---|---|

p1 | (x1, y1) | a(p1), b(p1), c(p1), d(p1) |

p2 | (x2, y2) | a(p2), b(p2), c(p2), d(p2) |

… | … | … |

Notice that each point has a different set of cubic coefficients and that these can be written as a long column vector. The above inline functions expect scalar (i.e. just one) values for `a`

, `b`

, `c`

, `d`

and will give an error if you tried to give it the `a`

, `b`

, `c`

and `d`

parameters as vectors.

To vectorize the expression, you use the `vectorize`

command:

R1 = inline(vectorize(char(r(1))),'a','b','c','d'); R2 = inline(vectorize(char(r(2))),'a','b','c','d'); R3 = inline(vectorize(char(r(3))),'a','b','c','d');

In this case, `vectorize`

is a simple string operation that puts a `.`

in front of every operator in the string that needs to be vectorized in order to work properly with vector `a`

, `b`

, `c`

and `d`

:

>> vectorize(char(r(1))) ans = 1./6./a.*(36.*c.*b.*a-108.*d.*a.^2-8.*b.^3+12.*3.^(1./2).*(4.*c.^3.*a-c.^2.*b.^2-18.*c.*b.*a.*d+27.*d.^2.*a.^2+4.*d.*b.^3).^(1./2).*a).^(1./3)- [...]

After this, all you need to do is set up your `a`

, `b`

, `c`

and `d`

as column vectors, with each row corresponding to values for each point you care about and then call `R1`

, `R2`

and `R3`

.

syms a b c d x r=solve(a*x^3+b*x^2+c*x+d); R1 = inline(vectorize(char(r(1))),'a','b','c','d'); R2 = inline(vectorize(char(r(2))),'a','b','c','d'); R3 = inline(vectorize(char(r(3))),'a','b','c','d'); % enter some values to play with. % note that there are 5 values for each coefficient. atemp = rand(5,1)*3; btemp = rand(5,1)*3; ctemp = rand(5,1)*3; dtemp = rand(5,1)*3; % each roots variable will have 5 elements, % which corresponds to each element in the coefficient % vectors. roots1 = R1(atemp,btemp,ctemp,dtemp); roots2 = R2(atemp,btemp,ctemp,dtemp); roots3 = R3(atemp,btemp,ctemp,dtemp);