Time to revive this blog. Compiling open cv 4.1 on a raspberry pi

Open CV 4.1 is out. This almost idiot proof tutorial should help you compile OpenCV 4 on a raspberry pi.

Step 1. Lets download the source for OpenCV 

https://opencv.org/releases/
I download version 4.1. If there's a newer version you can use that too

Step 2. Bring the Raspberry Pi up to par.

I mostly followed this blog https://towardsdatascience.com/setting-up-raspberry-pi-for-computer-vision-installing-opencv-e0f973d38819?gi=57ee64265a76

Skip the part about downloading OpenCV and just get the source file from step 1.

Save the following code to a shell script and do chmod +x your_script.sh 

In addition to the dependencies below to view images
sudo apt-get install libgtk2.0-dev pkg-config  

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

sudo apt-get update && sudo apt-get upgrade sudo apt-get install build-essential cmake pkg-config sudo apt-get install libjpeg-dev libtiff5-dev libjasper-dev libpng12-dev sudo apt-get install libavcodec-dev libavformat-dev libswscale-dev libv4l-dev sudo apt-get install libxvidcore-dev libx264-dev sudo apt-get install libgtk2.0-dev libgtk-3-dev sudo apt-get install libatlas-base-dev gfortran sudo apt-get -y remove x264 libx264-dev ## Install dependencies sudo apt-get -y install build-essential checkinstall cmake pkg-config yasm sudo apt-get -y install git gfortran sudo apt-get -y install libjpeg8-dev libjasper-dev libpng12-dev sudo apt-get -y install libtiff5-dev sudo apt-get -y install libtiff-dev sudo apt-get -y install libavcodec-dev libavformat-dev libswscale-dev libdc1394-22-dev sudo apt-get -y install libxine2-dev libv4l-dev cd /usr/include/linux sudo ln -s -f ../libv4l1-videodev.h videodev.h cd $cwd sudo apt-get -y install libgstreamer0.10-dev libgstreamer-plugins-base0.10-dev sudo apt-get -y install libgtk2.0-dev libtbb-dev qt5-default sudo apt-get -y install libatlas-base-dev sudo apt-get -y install libmp3lame-dev libtheora-dev sudo apt-get -y install libvorbis-dev libxvidcore-dev libx264-dev sudo apt-get -y install libopencore-amrnb-dev libopencore-amrwb-dev sudo apt-get -y install libavresample-dev sudo apt-get -y install x264 v4l-utils # Optional dependencies sudo apt-get -y install libprotobuf-dev protobuf-compiler sudo apt-get -y install libgoogle-glog-dev libgflags-dev sudo apt-get -y install libgphoto2-dev libeigen3-dev libhdf5-dev doxygen

Step 3. Lets build OpenCV using CMAKE

Before we do that make sure you do a pip3 install numpy
and sudo apt-get -y install cmake

Now you should be good to go. Create a build folder inside of open-cv

Inside this folder, make a sh file called build_me.sh 

Put the following code in the file and make it executable (chmod +x build_me.sh). I added something to include GTK, not sure if this will work but it's worth a try.

1	cmake ../ -D CMAKE_BUILD_TYPE=RELEASE \ -D CMAKE_INSTALL_PREFIX=/usr/local \ -D WITH_GTK=ON \ -D INSTALL_PYTHON_EXAMPLES=ON \ -D OPENCV_EXTRA_MODULES_PATH=~/opencv_contrib-4.1.0/modules \ -D BUILD_EXAMPLES=ON ..

Execute this file ./build_me.sh. It should make the installation files required for open-cv

Step 4. Lets make open-cv 

Before running make. Be sure to expand your swap space https://wpitchoune.net/tricks/raspberry_pi3_increase_swap_size.html

Do not run parallel make if you do not have a heat sink. It will over heat and freeze and you'll probably have to start the make over again.

Navigate to the cd build folder and run make

After make has completed, run sudo make install 

This builds all the install files you need

cd into python_loader and do a sudo python3 setup.py install 

Step 5. Lets check if it works

Bring up python in new terminal and do import cv2 if this works then you have open cv successfully installed.

Source panel method, solving potential flow around a body

Introduction

The source-panel is one of the simplest ways to model pressure along a body subjected to a flow at 0 angle of attack. The method ensures that flow cannot pass through the body - the normal component of the velocity on each panel is 0.

Source panel is not useful by itself for calculating lift

In order to use Potential Flow the following assumptions must be satisfied

- Incompressible

- Steady flow

- Low speed
- Irrotational

Approach

1. Introduction to Sources

2. Source sheet along a panel

3. Integration of a source sheet to calculate source strength

4. Effects of one panel's source sheet on the other and how to program the source-panel method

Sources/Sinks

A type of flowfield where the velocity is purely radial. You can think of it as a garden hose attached to a pan, the hose is turned on and fluid flow is coming out of the hose into the pan. The velocity potential $\phi$ is represented as $$ \phi = \frac{\gamma}{2\pi} log(r)$$

"r" is the distance to an arbitrary point by a source $\sqrt{(x-x0)^2 + (y-y0)^2}$

 The streamline $\psi$ is represented as $$\psi = \frac{\gamma}{2\pi} \theta$$ 

Velocity components u and v can both be derived from either the streamline or the velocity potential

$$u = \frac{\partial \phi}{\partial x} = \frac{\partial \psi}{\partial y}$$

$$v = \frac{\partial \phi}{\partial y} = -\frac{\partial \psi}{\partial x}$$

$$u = \frac{\gamma}{2\pi}\frac{x-x_0}{(x-x_0)^2 + (y-y_0)^2}$$

$$v = \frac{\gamma}{2\pi}\frac{y-y_0}{(x-x_0)^2 + (y-y_0)^2}$$

Panel Theory
The velocity potential around a body can be decomposed into $\phi = \phi_{\inf} + \phi_s + \phi_v$ freestream + source + vortex. For this example we will just consider the source and freestream.

$$ \phi_s = \int \frac{\gamma(s)}{2\pi} ln(r) ds $$ ds is the panel length and $q(s)$ is the source strength of the panel.

We don't immediately know $\gamma(s)$ this will be solved for in our system of equations.

$\gamma(s)$ depends on the free stream and the effects of all the other sources. The "r" term is the distance from panel i to panel j, j+2, j+3.

"r" is defined as $r_{i,j} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$, the radius from panel i to panel j.

In the Hess smith method, we assume a constant source throughout the panel. So in the figure above panel j's source directly effects panel i. One of the boundary conditions that we will assume is that there is no flow going through the center of panel i $V_i = 0$. We have to integrate the effects of panel j through the center of panel i.

$$
\phi_{i,j} = \int_j \frac{\lambda j}{2\pi} Ln(\sqrt{(x_i-x_j(s_j))^2+(y_i-y_j(s_j))^2})ds
$$

We want the boundary condition on panel i to be such that the velocity normal to the surface is 0. How do we achieve this? Take a look at the velocity potential. The velocity with respect to x is $\frac{\partial \phi}{\partial x}$, the velocity normal to the surface is similarly found to be $\frac{\partial \phi}{\partial n_i}$.

Finding $\frac{\partial \phi}{\partial n_i}$ seems challenging and it is, but we can break it down into the following.
$$ \frac{\partial}{\partial n_i} Ln(r_{i,j}) = \frac{1}{r_{i,j}} \frac{\partial r_{i,j}}{\partial n_i}$$
By chain rule:
$$
\frac{1}{r_{i,j}} \frac{\partial r_{i,j}}{\partial n_i} = \frac{1}{2 r_{i,j}} \frac{2(x_{ci}-x_j(s)) \frac{dx_i}{dn_i} + 2(y_{ci}-y_j(s)) \frac{dy_i}{dn_i}}{\sqrt{(x_{ci}-x_j(s))^2+(y_{ci}-y_{j}(s))^2}}
$$

Alternative View
It gets complicated but it is solvable. However, for me it's easier to view the x and y coordinates of panel j as a line, instead of using "s" as my parameter, I'll switch to using "t", t varies from 0 to 1

\begin{align}
s = \sqrt{(x_0+c_1 t)^2 + (y_0+v_t)^2} \\
ds = \frac{1}{2} \frac{2u(x_0 + c_1 t) + 2v(y_0 + c_2 t)}{\sqrt{(x_c-x(t))^2 + (y_c - y(t))^2}}
\end{align}

But it's still complicated. We can use the fact that $V \cdot \hat{n}$ is nothing but $ u\cdot \hat{n} + v \cdot \hat{n}$ Boundary conditions at "i" are $u \cdot \hat{n} = 0$ and $v \cdot \hat{n} = 0$ at $x_c$ and $y_c$. Note: $x_c$ and $y_c$ and $x_ci$ and $y_ci$ represent the location at the midpoint of panel i.

$$
\frac{1}{r_{i,j}} \frac{\partial r_{i,j}}{\partial n_i} = \frac{1}{2 r_{i,j}} \frac{2(x_{ci}-x_j(t)) \frac{dx_i}{dn_i} + 2(y_{ci}-y_j(t)) \frac{dy_i}{dn_i}}{\sqrt{(x_{ci}-x_j(t))^2+(y_{ci}-y_{j}(t))^2}}
$$

\begin{align}
x(t) = x_0 + c_1 t \\
y(t) = y_0 + c_2 t
\end{align}

From computer graphics we can find parameters $c_1$ and $c_2$ easily without using any angles.
\begin{align}
c_1 = x_1-x_0\\
c_2 = y_1-y_0
\end{align}
Now going back to the velocity potential, we have to integrate the velocity potential at "i" caused by all the sources in panel j from 1 to n where $j \neq i$ for all j. We have to sum up the effects of the source terms of all the panels!

$$
\Sigma_{j\neq i}^n V_{j}\cdot n_i= \Sigma_{j\neq i}^n \frac{\gamma_j}{2\pi} \int_{0}^{1} \frac{(x_{ci}-x(t)_j)cos(\beta_i) + (y_{ci}-y(t)_j)sin(\beta_i)}{(x_c-x(t)_j)^2 + (y_c - y(t)_j)^2} dt
$$

What do we do at i=j? Remember that with sources we have half in and half out.

We use the above information to help build a matrix Ax = b

$$A_{i,i} = 0.5$$

$$ A_{i,j} = \Sigma_{j\neq i} \frac{1}{2\pi} \int_{0}^{1} \frac{(x_{ci}-x(t)_j)cos(\beta_i) + (y_{ci}-y(t)_j)sin(\beta_i)}{(x_c-x(t)_j)^2 + (y_c - y(t)_j)^2} dt $$

We enforce the flow tangency condition ${\vec{V_{\inf}}} \cdot \hat{n_i}$ where $\vec{V_{\inf}} = [V_{\inf}cos(\alpha), V_{\inf}sin(\alpha)$, $n_i = [cos(\beta_i), sin(\beta_i)]$. By dot product and moving terms to rhs: 

$$ b_{i} =-\vec{V_{\inf}}cos(\alpha - \beta) $$

code coming soon

Transfinite Interpolation - Algebraic way of transposing computational domain to physical domain

Transfinite Interpolation was first described by William Gordon in 1973, it's used for transcribing the computational domain to the physical domain. I've been reading Handbook of Grid Generation by J.Thompson. Transfinite interpolation can be really confusing to people who are new to CFD and Grid work.

There is limited resources online on how Transfinite Interpolation really works. The goal of this post is to explain TFI in detail:

General Concept

From J.Thompson. Shows computational domain transformed into the physical domain

TFI comes from the Boolean Sum Formulation - J.Thompson Handbook of Grid Generation
\begin{align}
U(\xi,\eta,\zeta) = \Sigma_{i=1}^L \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \frac{\partial^n X(\xi_i, \eta, \zeta)}{\partial \xi^n}\\
V(\xi,\eta,\zeta) = \Sigma_{j=1}^M \Sigma_{m=0}^Q \beta_{j}^{m}(\eta) \frac{\partial^m X(\xi, \eta_j, \zeta)}{\partial \eta^m}\\
W(\xi,\eta,\zeta) = \Sigma_{k=1}^N \Sigma_{l=0}^R \gamma_{k}^{l}(\zeta) \frac{\partial^l X(\xi, \eta, \zeta_k)}{\partial \zeta^l}\\
\end{align}

$$X(\xi,\eta,\zeta)  =
\begin{bmatrix}
x(\xi,\eta,\zeta) \\
y(\xi,\eta,\zeta) \\
z(\xi,\eta,\zeta) \\
\end{bmatrix}$$

Tensor Products - not defined correctly in Thompson, but defined correctly here.
\begin{align}
UV = VU = \Sigma_{i=1}^L \Sigma_{j=1}^M \Sigma_{m=0}^Q \Sigma_{n=0}^R \alpha_{i}^{n}(\xi) \beta_{j}^{m}(\eta) \frac{\partial^{nm} X(\xi_i,\eta_j,\zeta)}{\partial \eta^{NM} \partial \xi^n}\\
UW = WU = \Sigma_{i=1}^L \Sigma_{k=1}^N \Sigma_{l=0}^R \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \gamma_{k}^{l}(\zeta) \frac{\partial^{ln} X(\xi_i,\eta,\zeta_k)}{\partial \zeta^{l} \partial \xi^{n}}\\
VW = WV = \Sigma_{j=1}^M \Sigma_{k=1}^N \Sigma_{m=0}^Q \Sigma_{l=0}^R \beta_{j}^{m}(\eta) \gamma_{k}^k (\zeta) \frac{\partial^{lm} X(\xi,\eta_j,\zeta_k)}{\partial \zeta^l \partial \eta^m}
UVW = \Sigma_{i=1}^L \ Sigma_{j=1}^M \Sigma_{k=1}^N \Sigma_{l=0}^R \ Sigma_{m=0}^Q \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \beta_{j}^{m}(\eta) \gamma_{k}^{l}(\zeta) \frac{\partial^{lmn} X(\xi_i,\eta_j,\zeta_k) }{\partial \zeta^{l} \partial \eta^m \xi^{n}}
\end{align}

$$ X(\xi_i,\eta_j,\zeta_k)  = U+V+W - UV-UW-VW-UVW$$

Looks confusing, lots of variables. The coefficients $\alpha$, $\beta$, and $\gamma$ depend on what method we use, they are called blending functions. Thompson goes into detail with them however, for the purpose of this post we will consider the simple case, linear interpolation.

Linear Interpolation the blending functions are simple, P=Q=R=0, L=M=N=2.
We are going to consider the 2D case
\begin{align}
\alpha_{1}^0 = 1-\xi \\
\alpha_{2}^0 = \xi \\
\beta_{1}^0 = 1-\eta \\
\beta_{2}^0 = \eta\\
U(\xi_i,\eta_j) = (1-\xi_i) X(0,\eta_j) + \xi_i X(1,\eta_j) \\
V(\xi_i,\eta_j) = (1-\eta_j)X(\xi_i,0) + \eta_j X(\xi_i,1)\\
UV(\xi_i,\eta_j) = (1-\xi_i)(1-eta_j) X(0,0) + \xi_i (1-\eta_j) X(1,0) + (1-\xi_i)(\eta_j) X(0,1) + \xi_i \eta_j X(1,1)
\end{align}

How is it used?
Take the simple duct example below. In order for us to map a computational domain onto the duct we have to convert this geometry from a function of x(x,y) y(x,y) to $x(\xi,\eta)$ and $y(\xi,\eta)$. Then we can use the above equations to map the computational domain.

In my code I do a simple mapping, nothing complicated. There may be a few bugs if you change geometries. $\xi$ and $\eta$ are passed into to a function and that outputs (x,y) coordinates. So where you see $X(\xi,\eta)$ the output is a vector $(x,y)$ then that vector is multiplied by $(1-\xi)$ or $\xi$ to get U,V, UV.
$$X(\xi,\eta) = U+V-UV$$

Grid mapped onto duct

Link to Matlab Code: GitHub

Finite Difference: 2nd Order Central, Any order

2nd order Central Difference can be found by adding the Forward and Backward Differences

Forward Difference:
\begin{equation}
u(x_0+\Delta x,y_0) = u(x_0, y_0) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} + \frac{\Delta x^3}{3!} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)
\end{equation}

Backward Difference:
\begin{equation}
u(x_0-\Delta x,y_0) = u(x_0, y_0) - \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} - \frac{\Delta x^3}{3!} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)
\end{equation}

Forward + Backward:
\begin{align}
u(x_0+\Delta x,y_0) + u(x_0-\Delta x,y_0) = 2u(x_0,y_0) + \frac{2 \Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} + \frac{2 \Delta x^4}{4!} \frac{\partial^4 u}{\partial x^4} \\
\frac{u(x_0+\Delta x,y_0) -2 u(x_0,y_0)+ u(x_0-\Delta x,y_0)}{\Delta x^2} = \frac{\partial^2 u}{\partial x^2} + O(\Delta x^2)
\end{align}

There's a pattern here. We can formulate finite difference schemes of any order given a set number of points. Take for example these 3 points. We can use these 3 points to find $u_x$ with $O(\Delta x^3)$ accuracy.

\begin{align} P(u_{i-1,j}) = u_{i,j} - \Delta x u_x + \frac{\Delta x^2}{2!} u_{xx} - \frac{\Delta x^3}{3!} u_{xxx} \\ Q(u_{i-2,j}) = u_{i,j} -2 \Delta x u_x + \frac{4 \Delta x^2}{2!} u_{xx} + \frac{3^3 \Delta x^3}{3!} u_{xxx} \end{align}

We want to set the term $u_{xx}=0$ and $u_x = 1$ we then have to solve the equation 

$$\frac{P}{2!} + \frac{4Q}{2!} = 0$$

$$P -2Q = 1$$

Substituting the coefficients $P=-2,Q=1/2$: 

$$-2 u_{i-1,j} + 2 u_{i,j} = -2\Delta x u_x + O(\Delta x^3)$$

$$1/2 u_{i-2,j} =  1/2 u_{i,j} - \Delta x u_x + O(\Delta x^3)$$

$$1/2 u_{i-2,j} - 2 u_{i-1,j} - 1.5 u_{i,j} = -3 \Delta x u_x + O(\Delta x^3)$$

We can use the same idea with 4 points and achieve a greater accuracy.

Finite Difference: Central Difference

The equation for Central Difference can be found by subtracting the equations for the Forward and Backward differencing schemes.

Forward Difference:
\begin{equation}
u(x_0+\Delta x,y_0) = u(x_0, y_0) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} + \frac{\Delta x^3}{3!} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)
\end{equation}

Backward Difference:
\begin{equation}
u(x_0-\Delta x,y_0) = u(x_0, y_0) - \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} - \frac{\Delta x^3}{3!} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)
\end{equation}

Subtracting these two equations $u(x_0+\Delta x,y_0) - u(x_0-\Delta x,y_0)$
\begin{align}
u(x_0+\Delta x,y_0) - u(x_0-\Delta x,y_0) = 0 + 2\Delta x \frac{\partial u}{\partial x} + 0 + O(\Delta x^3)\\
\frac{u(x_0+\Delta x,y_0) - u(x_0-\Delta x,y_0)}{2 \Delta x} = \frac{\partial u}{\partial x} + O(\Delta x^2)
\end{align}

Finite Difference: Forward and Backward Difference

When I first saw finite difference equation for the first derivative of a function "u" it seemed easy. It was just the derivative $\frac{\partial u}{\partial x} = \frac{u(x+\Delta x)}{\Delta x}$, but there's some form of art to how finite differencing really works.

Forward Difference:
\begin{equation}
\frac{\partial u}{\partial x} =\frac{u(x_0 + \Delta x, y_0) - u(x_0,y_0)}{\Delta x}
\end{equation}

Backward Difference:
\begin{equation}
\frac{\partial u}{\partial x} =\frac{ u(x_0,y_0) - u(x_0 - \Delta x, y_0)}{\Delta x}
\end{equation}

Two common schemes, we have the forward difference scheme and backward difference, but how did we get to here? In order to really understand finite difference we need to look at the Taylor Series.

\begin{equation}
u(x_0+\Delta x,y_0) = u(x_0, y_0) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} +...+ \frac{\Delta x^n}{n!} \frac{\partial^n u}{\partial x^n}
\end{equation}

Using the Taylor series we can derive the equation of the Forward Difference
\begin{equation}
\frac{u(x_0 + \Delta x, y_0) - u(x_0,y_0)}{\Delta x} = \frac{\partial u}{\partial x} + O(\Delta x)
\end{equation}

On the right hand side we have $O(\Delta x)$ this comes from the division on both sides by $\Delta x$; $\frac{O(\Delta x^2)}{\Delta x} = O(\Delta x)$ this is important when you're trying to find the order of accuracy of your finite difference equation.

To get backward differencing all we have to do is replace $\Delta x$ with $-\Delta x$

Poisson Equation Setup

The objective of this tutorial is to help you think of how to effectively program a solution to the Poisson Equation. The tutorial will go into how to think of node linking within matrices which can often be confusing.

Poisson Equation is defined to be $-\nabla{f} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}$ . Using central differencing the equation can be discretized to the following:
\begin{equation}
\frac{1}{h^2} \big( u(x_{i-1},y_{j})+u(x_{i},y_{j-1})-4u(x_{i},y_{j})+
u(x_{i+1},y_{j} +u(x_{i},y_{j+1}) \big) = f(x_{i}, y_{j})
\end{equation}

The grid is setup to be $5x5$, $n=4$. the formula for creating the grid should always be $(n+1)x(n+1)$. The total number of nodes will be $(n+1)^2$. We choose $(n+1)$ because it allows us to model the boundary.

Node list:

1	6	11	16	21
2	7	12	17	22
3	8	13	18	23
4	9	14	19	24
5	10	15	20	25

In the figure below Node 1 is connected to Node 6

We see that there is a pattern for node connectivity, N1 is connected to N6 and N2. N7 is connected to 2 and 12.

There is an easy way of programming this kind of connectivity. Take a look at the node table again, we see that if we iterate along the row by $i+1$ starting at node 1 we get 2. If we do the same thing except we iterate across the column by $j+1$ we arrive at 6. The same is true for node 7.

 function Poisson(n)  
 A = zeromatrix((n + 1)2; (n + 1)2)  
 F = zeros((n + 1)2,1)  
 G = reshape(1:((n + 1)2),n+1,n+1); . Note: 1:((n + 1)2) = 1,2,3,...,(n + 1)2 Reshape changes the  
 size of the vector 1,2,3,...,(n + 1)2 to a matrix n+1 by n+1  
 for i  0 to n do  
    for j  0 to n do  
   row = G(i+1,j+1);  
   if i=0 or j=0 or i=n or j =n then . We are on the boundary  
    A(row,row) = -4;  
    F(row,1) = 0;  
   else  
    A(row,G(i+1,j)) = 1;  
    A(row,G(i+1,j+2)) = 1;  
    A(row,G(i,j+1)) = 1;  
    A(row,G(i+2,j+1)) = 1;  
    A(row,G(i+1,j+1)) = -4;  
    F(row,1) = 1;  
 return A,F . Now solve for Ax=F  

120^2 x 120^2 Poisson Solution using Sparse Matricies