Vectorize symbolic expressions in MATLAB

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.

Problem Statement

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):

Find the nearest point on the parabola to the point in question

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)

Solving in MATLAB

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:

PointCoordinatesCubic 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.

Vectorize the Expression

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.

Complete Code Snippet

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);

Discussion

I would love to hear your feedback. Enter your comment below [ Terms of Use ]:
CTMES
 

About Peter Yu I am a research and development professional with expertise in the areas of image processing, remote sensing and computer vision. I received BASc and MASc degrees in Systems Design Engineering at the University of Waterloo. My working experience covers industries ranging from district energy to medical imaging to cinematic visual effects. I like to dabble in 3D artwork, I enjoy cycling recreationally and I am interested in sustainable technology. More about me...

Feel free to contact me with any questions about this site at [user]@[host] where [user]=web and [host]=peteryu.ca

Copyright © 1997 - 2019 Peter Yu