## Least squares fit of a surface to a 3D cloud of points in Python (with ridiculous application)

June 15, 2009 | categories: python, mathematics, estimation | View Comments

The floor in the room above the kitchen in our house has a floor that slopes almost 1 inch per foot for half of the room. Experimentally, we have found that this is steep enough to make a desk chair roll– kind of irritating, particularly to my special ladyfriend who happens to occupy such a chair in that zone.

To compensate for the slope, I decided to fit a stack of masonite sheets to the curve of the floor. Unfortunately, the floor slopes nonlinearly in two directions, like the rounded corner of a swimming pool. After making a series of measurements of the floor, I decided to fit a polynomial in two variables to the cloud of points using a least squares estimate.

(In retrospect, the floor was close enough to singly-curved that I could have gotten away with a linear fit.) The blue thing in the picture is a level. I shimmed the level until it was worthy of its name, and then measured the distance to the floor with calipers.

**Estimating the flatness of the floor**

I’ll recount the basics from my earlier post about least squares fitting in Python. Skip ahead to the next section if you read that already.

As before, the first step is to arrange the equations in canonical form: Ax=y where:

```
* A is an M x N full rank matrix with M > N (skinny)
* x is a vector of length N
* y is a vector of length M
```

As matrices, this looks like

In polynomial fitting, A is called the Vandermonde matrix and takes the form:

**The 3D case**

In the 2D case, we’re trying to find polynomial in x such that f(x) approximates y. In the 3D case at hand, we have two independent variables, so we’re looking for a polynomial in x and y such that f(x, y) approximates z. Rather than the 2D case:

we want the final output to look like this:

(Spike Curtis astutely notes in the comments that I am omitting the cross terms, such as xy, and refers us to a Matlab script. Spike is right, but the floor’s already fixed.)

We can use two Vandermonde matrices next to each other.

Here’s the Python code that creates the two Vandermonde matrices and joins them into one matrix. x, y, and z are lists of corresponding coordinates, so, for example, x[5], y[5] and z[5] are the coordinates of one point that the surface should approximate. The order of the points is not important.

import numpy as np z = [0.0, 0.695, 1.345, 1.865, 2.225, 2.590, 0.0, 0.719, 1.405, 1.978, 2.398, 2.730, 0.0, 0.789, 1.474, 2.064, 2.472, 2.775, 0.0, 0.763, 1.453, 1.968, 2.356, 2.649] x = [0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0] y = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 48.0, 48.0, 48.0, 48.0, 48.0, 48.0, 72.0, 72.0, 72.0, 72.0, 72.0, 72.0] degree = 3 thickness = 0.167 # Set up the canonical least squares form Ax = np.vander(x, degree) Ay = np.vander(y, degree) A = np.hstack((Ax, Ay)) # Solve for a least squares estimate (coeffs, residuals, rank, sing_vals) = np.linalg.lstsq(A, z) # Extract coefficients and create polynomials in x and y xcoeffs = coeffs[0:degree] ycoeffs = coeffs[degree:2 * degree] fx = np.poly1d(xcoeffs) fy = np.poly1d(ycoeffs)

Once I knew the coefficients of the polynomial approximation, I could calculate the contours of the masonite. From a measurement of the stack of masonite, I knew the average thickness is 0.167 inches. (Thanks to Dr. Alex T. Tung, Ph.D., for helping me get the masonite home from Home Depot.) To find out where the first layer should end, I picked a series of stations spaced every 12 inches along the length of the floor in the y-direction. Along those stations, I solved f(x, y) = 0.167, f(x, y) = 2 * 0.167, f(x, y) = 3 * 0.167 and so forth.

In practice, solving f(x, y) = c, where c is a constant, means finding the roots of the equation f(x, y) - c = 0. (The mathematicians call this solving the homogeneous equation.) In Python, the numpy.roots method solves the homogeneous case. For each contour/section crossing, I generated a polynomial of the form f(x, y) - c and solved it with numpy.roots.

ystations = range(0, 84, 12) sections = [[np.poly1d(xcoeffs - [0,0,zoffset - fy(ypos)]) for zoffset in np.arange(thickness, max(z), thickness).tolist()] for ypos in ystations] pts = [[min(func.roots) for func in list_of_fs] for list_of_fs in sections]

For fabrication, I printed out a list of the locations where the masonite contours crossed the stations.

for (pt_list, ystation) in zip(pts, ystations): print('\nBoundaries at station y = {0} inches:'.format(ystation)) print('\t'.join(['{0:.3}'.format(pt) for pt in pt_list]))

Armed with my list of measurements, I headed to the garage and set up some sawhorses with a sheet of plywood to keep the masonite from bowing and flopping around. It took a few hours of marking points and cutting gentle curves with a jigsaw, but the results were delightful.

*The masonite platform (non-impressive view)*

*Level and gleeful*

*Side view of the masonite platform*

*Another side view of the masonite platform*

As you can see above, I haven’t fastened the layers together yet. They seem to be sticking together reasonably well.

In the end, I think the slope went from about a 2.75″ drop across the 48″ span, to within 1/8″ of flat across the same distance. Sharon was delighted, and I found the whole project both more time-consuming and more satisfying than I expected.

## Calculating solar panel shading in Python

March 16, 2009 | categories: python, mathematics, solar | View Comments

Eventually, I'd like to install solar panels on our house, but I want to know that it will be worth it before I commit the money. Back in 2007, I wrote a library for calculating the path of the sun given the time and your location on the earth. Since then, I've been thinking about the next step, which is calculating how much neighboring houses and trees would obstruct the sun at different times of the day. In a nutshell, I wanted to replicate the calculations performed by devices like the Solmetric Suneye without paying $1500 for the device. To do the job right, I still need a quality circular fisheye lens ($680 new), but I've got a first approximation working with a $13 door peephole from McMaster.

(I do realize that I could just get a solar contractor to do all the site evaluation for me for free. That misses the point. If all I wanted was some sort of fiscal efficiency, I'd get a job as a financial executive or maybe something different from that, like a bank robber.)

Back to the task at hand.

Below is the image I started with, which I took back in the fall, before the leaves fell off the trees. This was using a Canon SD450 camera and the aforementioned peephole. I didn't take particular care to level the camera or center the peephole.

Using the Python Imaging Library, I cropped the photo to be square and switched to black and white mode. Then I mapped concentric rings in the image below to the subsequent straightened image.

Here's the straightening code. The guts are in the nested for loops near the end. (If there are any Pythonistas out there who know how to iterate over concentric rings using a list comprehension or map(), please let me know. The code below is functional and reasonable clear, but a little slow.)

# Code under GPLv3; see pysolar.org for complete version. def despherifyImage(im): (width, height) = im.size half_width = im.size[0]/2 half_height = im.size[1]/2 inpix = im.load() out = Image.new("L", (width, half_height)) outpix = out.load() full_circle = 1000.0 * 2 * pi for r in range(half_width): for theta in range(int(full_circle)): (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width) (outx, outy) = (width - width * (theta/full_circle) - 1, r) outpix[outx, outy] = inpix[inx, iny] return out

The straightening works pretty well, but there is a little distortion that I think is caused by the peephole not being a perfectly spherical lens. Also, because the peephole is not centered on the camera lens, my concentric ring transformation isn't centered either.

From here, I needed to figure out where the sky ends and the buildings or trees start. Unfortunately, the buildings and trees are both lighter *and* darker than the sky in different places, so I can't just look for one type of transition. To detect the edge, I scan down each column of pixels and calculate the difference in darkness between consecutive pixels, ignoring whether the change was from light to dark or the reverse. The image below shows those differences, amplified by 10x to increase the contrast. The mathematicians call this calculation a finite forward difference.

Here's the finite difference code.

# Code under GPLv3; see pysolar.org for complete version. def differentiateImageColumns(im): (width, height) = im.size pix = im.load() for x in range(width): for y in range(height - 1): pix[x, y] = min(10 * abs(pix[x, y] - pix[x, y + 1]), 255) return im

The last step is to scan down each column looking for the first large value. The first change that crosses a threshold is recorded. The final output is an array of values that measure the angle of the highest obstruction as a function of direction. As a sanity check, I drop a red dot on each value. It's hard to make out in the thumbnail below, but if you click on the image below, you'll get a larger version where you can see the red dots work pretty well.

I think the next step will be to calculate the total energy delivered per year using Pysolar.

Here's the full code.

# Code under GPLv3; see pysolar.org for complete version. from PIL import Image from math import * import numpy as npdef squareImage(im): (width, height) = im.size box = ((width - height)/2, 0, (width + height)/2, height) return im.crop(box)

def despherifyImage(im): (width, height) = im.size half_width = im.size[0]/2 half_height = im.size[1]/2 inpix = im.load() out = Image.new("L", (width, half_height)) outpix = out.load() full_circle = 1000.0 * 2 * pi for r in range(half_width): for theta in range(int(full_circle)): (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width) (outx, outy) = (width - width * (theta/full_circle) - 1, r) outpix[outx, outy] = inpix[inx, iny] return out

def differentiateImageColumns(im): (width, height) = im.size pix = im.load() for x in range(width): for y in range(height - 1): pix[x, y] = min(10 * abs(pix[x, y] - pix[x, y + 1]), 255) return im

def redlineImage(im): (width, height) = im.size pix = im.load() threshold = 300 for x in range(width): for y in range(height - 1): (R, G, B) = pix[x, y] if R + G + B > threshold: pix[x, y] = (255, 0, 0) break return im

im = Image.open('spherical.jpg').convert("L") im = squareImage(im)

lin = despherifyImage(im) d = differentiateImageColumns(lin).convert("RGB") r = redlineImage(d)

r.show()

## Least squares polynomial fitting in Python

January 24, 2009 | categories: python, mathematics, estimation | View Comments

A few weeks ago at work, in the course of reverse engineering a dissolved oxygen sensor calibration routine, Jon needed to fit a curve to measured data so he could calculate calibration values from sensor readings, or something like that. I don’t know how he did the fitting, but it got me thinking about least squares and mathematical modeling.

The tools I’ve used for modeling (Matlab, Excel and the like) and the programming languages I’ve used for embedded systems (C and Python, mostly) are disjoint sets. Who wants to use Matlab for network programming, when you could use Python instead? Who wants to use C for calculating eigenvalues, when you could use Matlab instead? Generally, nobody.

Unfortunately, sometimes you need to do both in the same system. I’m not sure that I’m making a wise investment, but I think that Python’s numerical libraries have matured enough that Python can straddle the divide.

So, here’s what I’ve figured out about least squares in Python. My apologies to most of you regular readers, as you will find this boring as hell. This one’s for me and the internet.

**Least squares estimation**

Generally, least squares estimation is useful when you have a set of unknown parameters in a mathematical model. You want to estimate parameters that minimize the errors between the model and your data. The canonical form is:

Ax=y where:

```
* A is an M \times N full rank matrix with M > N (skinny)
* x is a vector of length N
* y is a vector of length M
```

As matrices, this looks like

You want to minimize the L2-norm, ||Ax-y||. To be explicit, this means that if you define the residuals r = Ax - y, you want to choose the that minimizes sqrt{r_1^2 + ... r_M^2}. This is where the “least squares” name comes from– you’re going to choose the that results in the squares of r being least.

Once you know A and y, actually finding the best x (assuming you have a computer at your disposal) is quick. The derivation sets the gradient of ||r||^2 with respect to x to zero and then solves for x. Mathematicians call the solution the Moore-Penrose pseudo-inverse. I call it “what the least squares function returned.” In any case, here it is:

Fortunately, one need not know the correct name and derivation of a thing for it to be useful.

So now, how do you go about using this technique to fit a polynomial to some data?

**Polynomial fitting**

In the case of polynomial fitting, let’s say you have a series of data points, (t_i, y_i) that you got from an experiment. The t_i are the times at which you made your observations. Each y_i is the reading from your sensor that you noted at the time with the same i.

Now you want to find a polynomial y = f(t) that approximates your data as closely as possible. The final output will look something like

where you have cleverly chosen the x_i for the best fit to your data. You’ll be able to plug in whatever value of t you want and get an estimate of y in return.

In the conventional form, we are trying to find an approximate solution to Ax=y

In this case, A is called the Vandermonde matrix; it takes the form:

(This presumes that you want to use powers of t as your basis functions– not a good idea for all modeling tasks, but this blog post is about polynomial fitting. Another popular choice is sinusoidal functions; perhaps Google can tell you more on that topic.)

(Also, the Vandermonde matrix is sometimes defined with the powers ascending from left to right, rather than descending. Descending works better with Python because of the coefficient order that numpy.poly1d() expects, but mathematically either way works.)

(Back to polynomial fitting.)

Each observation in your experiment corresponds to a row of A; A_i * x = y_i. Presumably, since every experiment has some noise in it, there is no x that fits every observation exactly. You know all the t_i, so you know A.

As a reminder, our model looks like this:

You can check that A makes sense by imagining the matrix multiplication of each row of A with x to equal the corresponding entry in y. For the second row of A, we have y_2 = x_Nt_2 + . . . + x_2t_2^2 + x_1t_2^1 + x_0

which is just the polynomial we’re looking to fit to the data.

Now, how can we actually do this in Python?

Python code

Let’s suppose that we have some data that looks like a noisy parabola, and we want to fit a polynomial of degree 5 to it. (I chose 5 randomly; it’s a stupid choice. More on selection of polynomial degree in another post.)

import numpy as np import scipy as sp import matplotlib.pyplot as pltdegree = 5

# generate a noisy parabola t = np.linspace(0,100,200) parabola = t2 noise = np.random.normal(0,300,200) y = parabola + noise

So now we have some fake data– an array of times and an array of noisy sensor readings. Next, we form the Vandermonde matrix and find an approximate solution for x. Note that numpy.linalg.lstsq() returns a tuple; we’re really only interested in the first element, which is the array of coefficients.

# form the Vandermonde matrix A = np.vander(t, degree) # find the x that minimizes the norm of Ax-y (coeffs, residuals, rank, sing_vals) = np.linalg.lstsq(A, y) # create a polynomial using coefficients f = np.poly1d(coeffs)

Three lines of code later, we have a solution!

I could have also used numpy.polyfit() in place of numpy.linalg.lstsq(), but then I wouldn’t be able to write a blog post about regularized least squares later.

Now we’ll plot the data and the fitted curve to see if the function actually works!

# for plot, estimate y for each observation time y_est = f(t) # create plot plt.plot(t, y, '.', label = 'original data', markersize=5) plt.plot(t, y_est, 'o-', label = 'estimate', markersize=1) plt.xlabel('time') plt.ylabel('sensor readings') plt.title('least squares fit of degree 5') plt.savefig('sample.png')

This post took two-thirds of eternity to write. Feel free to ask questions or point out how confused I am below.

Next: regularized least squares! Or maybe selection of polynomial degree! Or maybe forgetting the blog and melting brain with TV!