Skip to content

Polynomial Fitting

import numpy as np
import matplotlib.pyplot as plt

Goal

Since an image is a discrete mapping of pixels, while the real world is "continuous", we want to fit the pixels by a continuous functions so that we can extract information.

Sliding Window Algorithm

  • Define a "pixel window" centered at pixel (x,y)(x, y) so that an 1D patch can be (0,y)(2x,y)(0, y)\sim (2x,y)
  • Fit a n-degree polynomial to the intensities (commonly n2n \leq 2)
  • Assign the poly's derivatives at x=0x=0 to pixel at windows's center
  • Slide window one pixel over until window reaches border

Taylor-series Approximation of Images

Consider a 1D patch from (0,y)(0,y) to (2x,y)(2x,y)
Consider the Taylor expansion of II centered at 00

I(x)=i=0diIdxixiI(0)i=0NdiIdxixiI(0)+RN+1(x)I(\vec x) = \sum_{i=0}^\infty \frac{d^iI}{d\vec x^i}\vec x^i I(0) \approx \sum_{i=0}^N \frac{d^iI}{d\vec x^i}\vec x^i I(0) + R_{N+1}(x)

Thus, we can write it as

I(x)=[1,x,x2/2,x3/6,...,xn/n!][I(0),dxI(0),dx22I(0),...,dxnnI(0)]TI(x) = [1, x, x^2/2, x^3/6,...,x^n/n!][I(0), d_xI(0), d^2_{x^2}I(0), ..., d^n_{x^n}I(0)]^T

Note that we have 2x+12x+1 pixels in the patch, ranges from x-x to xx, let their intensities be IxI_x. So that we have

[IxIx+1...Ix1Ix]=[1x...(x)n/n!1(x+1)...(x+1)n/n!............1(x1)...(x1)n/n!1x...xn/n!][I(0)dxI(0)...dxn1n1I(0)dxnnI(0)] \begin{bmatrix} I_{-x}\\ I_{-x+1}\\ ...\\ I_{x-1}\\ I_{x} \end{bmatrix}= \begin{bmatrix} 1 &-x &...&(-x)^n/n!\\ 1 &(-x+1) &...&(-x+1)^n/n!\\ ... &... &...&...\\ 1 &(x-1) &...&(x-1)^n/n!\\ 1 &x &...&x^n/n!\\ \end{bmatrix} \begin{bmatrix} I(0)\\d_xI(0)\\...\\d^{n-1}_{x^{n-1}}I(0)\\d^n_{x^n}I(0)\end{bmatrix}
I(2x+1)×1=X(2x+1)×ndn×1I_{(2x+1)\times 1}=X_{(2x+1)\times n}d_{n\times 1}

where I,XI, X are known and we want to solve dd

Since this is not always have a solution. We want to minimize the "fit error" IXd2=nIi2\|I-Xd\|^2 = \sqrt{\sum^n I_i^2}

When n=0n = 0, obviously dd is minimized at mean, i.e. 2x+1Ii/(2x+1)\sum^{2x+1} I_i / (2x+1)

Weighted least square estimation

If we want to extract more information, or in favor of, the neighboring points, we can have I=ΩII' = \Omega I where Ω\Omega is the weight function.

In 1-D case, let Ω=diag(ω1,...,ω2x+1)\Omega = diag(\omega_1, ..., \omega_{2x+1}), we want to solve ΩI=ΩXd\Omega I = \Omega X d which is to minimize Ω(IXd)2\|\Omega(I-Xd)\|^2

Random Sample Consensus (RANSAC)

Suppose we have some outliers, then the line fitting with least square is not a good solution.

x = np.arange(-10, 10, 0.5)
y = x * 0.5 + np.random.normal(0, 1, x.shape)
y += (np.random.randn(y.shape[0]) < -1.5) * 8
y -= (np.random.randn(y.shape[0]) < -1.5) * 8
fit = np.polyfit(x, y, 1)
f = fit[0] + fit[1] * x
plt.scatter(x, y)
plt.plot(x, f);

png

Algorithm

Given:
n=n= degree of freedom (unkonwn polynomial coefficients)
p=p= fraction of inliers
$t = $ fit threshold
$P_s = $ success probability
K=K= max number of iterations

RANSAC
for _ in K:
    randomly choose n + 1 pixels
    fit a n-degree polynomial
    count pixels whose vertical distance from poly is < t
    if there are >= (2w+1) * p pixels, mark them as inliers:
        do n-degree poly fitting on inliers 
        return the fitting
# if there is insufficient inliers
return None

Consider the value for KK:
Suppose inliers are independent distributed, hence

P(choose n+1 inliers)=pn+1P(\text{choose n+1 inliers}) = p^{n+1}
P(at least 1 outliers)=1pn+1P(\text{at least 1 outliers}) = 1 - p^{n+1}
P(at least 1 outliers in all K trials)=(1pn+1)KP(\text{at least 1 outliers in all K trials}) = (1-p^{n+1})^K
Ps=1(1pn+1)KK=log(1Ps)log(1pn+1)P_s = 1 - (1-p^{n+1})^K\Rightarrow K = \frac{\log(1-P_s)}{\log(1-p^{n+1})}

2D Image Patch

Given an image patch I(1:W,1:H),I[0,1]I(1:W, 1:H), I\in[0,1] is the intensity mapping

2D Taylor Series Expansion

Let I(0,0)I(0,0) be the center

I(x,y)=I(0,0)+xxI(0,0)+yyI(0,0)+12(x22Ix2(0,0)+y22Iy2(0,0)+xy2Ixy(0,0))+...\begin{align*} I(x,y) &= I(0,0) + x\partial_xI(0,0) + y\partial_yI(0,0) \\ & + \frac{1}{2}(x^2\frac{\partial^2I}{\partial x^2}(0,0) + y^2\frac{\partial^2I}{\partial y^2}(0,0) + xy\frac{\partial^2I}{\partial x\partial y}(0,0)) + ... \end{align*}

Directional Derivative

Using 1st order Taylor series approximation.
Let θ\theta be the direction, v=(cosθ,sinθ)v=(\cos \theta, \sin\theta) be the unit vector in direction θ\theta.
Consider the points (tcosθ,tsinθ)(t\cos\theta, t\sin\theta), i.e. tt units away from (0,0)(0,0) in the direction of θ\theta

I(tcosθ,tsinθ)=I(0,0)+I(0,0)+tcosθxI(0,0)+tsinθyI(0,0)I(t\cos\theta, t\sin\theta) = I(0,0) + I(0,0) + t\cos\theta\partial_xI(0,0) + t\sin\theta\partial_yI(0,0)

Then consider the directional derivative

vI(0,0)=limt0I(tcosθ,tsinθ)I(0,0)tcosθxI(0,0)+sinθyI(0,0)=[xI(0,0)yI(0,0)][cosθsinθ]\begin{align*}\partial_v I(0,0) &= \frac{\lim_{t\rightarrow 0}I(t\cos\theta, t\sin\theta) - I(0,0)}{t}\\ &\approx \cos\theta\partial_xI(0,0) + \sin\theta\partial_yI(0,0)\\ &=\begin{bmatrix}\partial_xI(0,0)&\partial_yI(0,0)\end{bmatrix}\begin{bmatrix}\cos\theta\\\sin\theta\end{bmatrix} \end{align*}

Let I=(xI,yI)\nabla I = (\partial_xI, \partial_yI). Note that

vI(0,0)={1v is parallel to I(0,0)0vI(0,0)(isophote)\partial_vI(0,0) = \begin{cases}1 &v\text{ is parallel to }\nabla I(0,0) \\0 &v\perp \nabla I(0,0) (\text{isophote})\end{cases}

Therefore, the gradient I\nabla I tells the direction and magnitude of max changes in intensity.

Discrete Derivative

Noting that the image are pixels, hence we can approximate the derivative by its limit definition.

xI(x,y)=12(I(x+1,y)f(x1,y))\partial_xI(x,y) = \frac{1}{2}(I(x+1,y) - f(x-1,y))
yI(x,y)=12(I(x,y+1)f(x,y1))\partial_yI(x,y) = \frac{1}{2}(I(x,y+1) - f(x,y-1))

Or we can use forward difference (x+1,x)(x+1, x) or backward difference (x,x1)(x,x-1) vs. (x1,x+1)/2(x-1,x+1)/2

Magnitude I(x,y)=xI(x,y)2+yI(x,y)2\|\nabla I(x,y)\| = \sqrt{\partial_xI(x,y)^2 + \partial_yI(x,y)^2} is the L2L_2 norm
Direction Θ(x,y)=arctan(xI/yI)\Theta(x,y) = \arctan(\partial_xI / \partial_yI)

Laplacian

Based on central difference, the approximation to the laplacian of II is

d2fdx2f(x+1)2f(x)+f(x1)\frac{d^2f}{dx^2} \approx f(x+1)-2f(x) + f(x-1)

For 2D functions,

2fx2f(x+1,y)2f(x,y)+f(x1,y)\frac{\partial^2f}{\partial x^2} \approx f(x+1, y)-2f(x,y) + f(x-1,y)

Hence

2f=2fx2+2fy2=f(x+1,y)2f(x,y)+f(x1,y)+f(x,y+1)2f(x,y)+f(x,y1)=f(x+1,y)+f(x1,y)+f(x,y+1)+f(x,y1)4f(x,y)=[010141010]\begin{align*} \nabla^2 f &= \frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}\\ &= f(x+1, y)-2f(x,y) + f(x-1,y) + f(x, y+1)-2f(x,y) + f(x,y-1)\\ &= f(x+1, y) + f(x-1,y) + f(x, y+1)+ f(x,y-1)-4f(x,y)\\ &= \begin{bmatrix} 0&1&0\\1&-4&1\\0&1&0 \end{bmatrix} \end{align*}