OpenCV feature detection and tracking

Wrote some utter shit code to get started with tracking via feature detection. Probably not the best way to go about the business.

Click on the screen to pick the green dot to follow

Uses the ORB feature finder. Each feature then has a descriptor and does brute force search to find best one.

Needs a lot more to be useful for anything. Suggestions: reject outliers, take best k matches and pick closest one to previous match, maintain a set of best current descriptors.

This is another approach to following something, but I’m not sure it is a great one.

import cv2
import numpy as np

orb = cv2.ORB()
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)

minpoint = []
mindes = []
def findNearestKeyPoint(x,y,kp,des):
    global minpoint, mindes
    mindist = 100000000000
    minpoint = kp[0]
    mindes = des[0]
    index = 0
    for point in kp:
        dist = ([0]-x)**2+([1]-y)**2
        if dist < mindist:
            mindist = dist
            minpoint = point
            mindes = des[index]
            #print des[index]
        index += 1

def service_mouse(event,x,y,flags,param):
    global kp, des
    if event == cv2.EVENT_LBUTTONDOWN:

cap = cv2.VideoCapture(0)



    _, frame =

    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    kp = orb.detect(gray,None)
    kp, des = orb.compute(gray, kp)
    img2 = cv2.drawKeypoints(frame,kp,color=(0,255,0), flags=0)
    if minpoint:, (int([0]),int([1])) , 4, [0,0,255],-1 )

        matches = bf.match(np.array([mindes]),des)
        # suggestions: get list of best, pick closest to knwon position, update descriptors

        matches = sorted(matches, key = lambda x:x.distance)
        newguy = kp[matches[0].trainIdx].pt, (int(newguy[0]),int(newguy[1])) , 10, [0,0,255],-1 )


    k = cv2.waitKey(5) & 0xFF
    if k == 27:





Knife Bot

Knife bot. Hide your children and wives and dogs.

It’s a hot glued arrangement of servos and home depot painter sticks (they’re free!)


We used python and numpy to solve for the needed angles as a function of the end position. We were having trouble with getting angles that are in the -np.pi/2 to np.pi/2 range until the suggestion of using [0,np.pi/2] as the initial position. In hindsight, a more robust method would be to make a loop attempting a set of reasonable starting positions, and accept them if they were attainable. Or to attempt a random set of initial conditions. Your choice. For only 2 variables, the computer will just chew that problem up. If you had 10 actuators, maybe you’d have more problems.

Turns out random angles mode was more fun though.

from scipy import optimize
import numpy as np
import serial

device = '/dev/cu.usbmodem1411'
L = 1.

def pos(theta, x):
    rho = L* (np.sin(theta[0]) + np.sin(theta[0] + theta[1]))
    z = L * (np.cos(theta[0]) + np.cos(theta[0] + theta[1]))
    return np.array([rho - x[0], z - x[1]])

def objective(theta,x):
    temp = pos(theta,x)
    return temp[0]**2+temp[1]**2

anglerange = (-np.pi/2,np.pi/2)
bnds = (anglerange,anglerange)

#cons = [{'type':'ineq', 'fun': con},
#        {'type':'eq', 'fun': con_real}]

def cleanAngle(angle):
    angle = (angle + np.pi) % (2 * np.pi) - np.pi
    return angle

def findAngles(x,y,z):

    rho = np.sqrt(x**2 + y**2)
    theta0 = np.arctan2(y,x)
    sol = optimize.root(lambda theta: pos(theta, [rho,z] ), [0,np.pi/2])
    #sol = optimize.minimize(lambda theta: objective(theta, [rho,z]), np.array([0,0]), method='SLSQP', bounds=bnds)
    print sol.success
    theta = sol.x
    theta = cleanAngle(theta)
    return [theta0, theta[0], theta[1]]

ser = serial.Serial(device, 19200, timeout=1)

import signal
import sys
def signal_handler(signal, frame):
        print('You pressed Ctrl+C!')
signal.signal(signal.SIGINT, signal_handler)

import time
while True:
    x = float(input('x: '))
    y = float(input('y: '))
    z = float(input('z: '))
    theta = findAngles(x,y,z)
    ser.write(str(theta[0])+ ',')
    ser.write(str(theta[1])+ ',')
    ser.write(str(theta[2])+ ';\r\n')
while True:
    theta = np.random.rand(3) * np.pi - np.pi/2
    ser.write(str(theta[0])+ ',')
    ser.write(str(theta[1])+ ',')
    ser.write(str(theta[2])+ ';\r\n')


This is the arduino code. Servos are on 9,10,11

We dialed in the Zero and Ration numbers by iteratively try to get angle 0 and 90 degrees manually over the serial monitor for each servo in turn.

The command strings is in the format


Comma separated angles for each motor in radians.



#include <Servo.h>
//#include <SoftwareSerial.h>
#define RATIO 425/(PI/4)
#define BRATIO -400/(PI/4)
#define CRATIO -550/(PI/4)
#define AZERO 1606
//#define AZERO 2020
#define BZERO 1645
#define CZERO 2560
#define D 200.

Servo servoA, servoB, servoC;  // create servo object to control a servo
String v = "";
int sign = 1;
float a = 0.;
float b = 0.;
float c = 0.;
float x_temp;
float y_temp;
float angles[3] = {0.,0.,0.};

//SoftwareSerial ser(2,3);

void setup()
  servoC.attach(11);  // attaches the servo on pin 9 to the servo object

void loop() {
  // put your main code here, to run repeatedly:
  if (Serial.available()) {
    char ch =;
    switch(ch) {
      case ';':
        for(int i = 0; i<3 ; i++){
        v += ch;

void resetV() {
  v = "";

void getAngles(String input, float *angles) {
  int comma_1 = 0;
  int comma_2 = 0;
  for (int i = 0; i < input.length(); i++) {
    if (input.substring(i, i+1) == ",") {
      if (comma_1 == 0) {
        comma_1 = i;
      } else {
        comma_2 = i;
  angles[0] = input.substring(0, comma_1).toFloat();
  angles[1] = input.substring(comma_1+1, comma_2).toFloat();
  angles[2] = input.substring(comma_2+1, input.length()).toFloat();

void move(float* angles) {

      servoA.writeMicroseconds(floor(angles[0]*RATIO) + AZERO);
      servoB.writeMicroseconds(floor(angles[1]*BRATIO) + BZERO);
      servoC.writeMicroseconds(floor(angles[2]*CRATIO) + CZERO);



More Opencv

Canny finds edges. Edges are basically places of very large derivative in the image. Then the canny algorithm cleans it up a little.

FindContours seems to be a clutch dude

This guy uses it to

import cv2
import numpy as np

cap = cv2.VideoCapture(0)


    # Take each frame
    _, frame =

    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    edges = cv2.Canny(gray,100,200)
    #ret,thresh = cv2.threshold(gray,50,255,cv2.THRESH_BINARY)

    kernel = np.ones((5,5),np.uint8)
    edges = cv2.dilate(edges,kernel,iterations = 5) # really chunks it up

    contours,hierarchy= cv2.findContours(edges,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)[-2:]
    cv2.drawContours(frame, contours, -1, (0,0,255), 3)
    k = cv2.waitKey(5) & 0xFF
    if k == 27:


The dilation reduces the number of contours to something more reasonable

Each contour is of the format [[[x y]], [[x y]], [[x y]]]

I had a problem with draw contours until I found out I needed to write onto a color image with it.


The good features to track

import cv2
import numpy as np

cap = cv2.VideoCapture(0)


    _, frame =

    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    corners = cv2.goodFeaturesToTrack(gray,25,0.01,10)
    corners = np.int0(corners)
    for i in corners:
        x,y = i.ravel(),(x,y),8,[0,0,255],-1) #image center radius color thickness

    k = cv2.waitKey(5) & 0xFF
    if k == 27:

The 25 is number of corners, 0.01 is a quality cutoff (1% of best corner quality found), 10 is minimum distance between corners


Ubuntu Virtualbox

Download an ubuntu install distro

Download Virtualbox

Make a new virtual machine. Linux 64-bit probably

Under storage go the the cd and click on it. Select the installation disk.

Decide to erase and install, it’s fine.

Insert the guest additions cd using the menu bar under devices. Install them. Now it should go full screen when you click the little green button on macs. You can slide all your fingers to switch workspaces.

TO be able to access shared folders

sudo adduser USERNAME vboxsf

Also, checking 3d acceleration under Settings > Display and giving it all the video ram you can seems to make it run much zippier.

I also gave it 3gb of ram because i had a little extra to spare. Seems to like it.

sudo apt-get install lxde

Can use a more lightweight desktop environment by selecting that little emblem at the login screen.

ipywidgets: useful little buggers

So, the matplotlib slider functionality is basically garbage.

Let’s see if ipywidgets and ipython notebooks do better

sudo pip install ipywidgets


ipython notebook


A basic program with sliders. You can go from there.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import interactive
from IPython.display import display

def beat_freq(f1=1.0, f2=1.0):
    times = np.linspace(0,3,100)
    signal = np.sin(2*np.pi*f1*times) + np.sin(2*np.pi*f2*times)
    #print(f1, f2, abs(f1-f2))
    #display(Audio(data=signal, rate=rate))

v = interactive(beat_freq, f1=(0.0,1.0), f2=(0.0,1.0))


Performing Some Laplace Experiments

So I’m looking at the sparse solver capabilities of python.

Trying a 3d poisson problem (electrostatics of a point charge) in V=0 boundary conditions

\nabla^2 \phi = \delta^3 (x)


from scipy import sparse
from scipy.sparse import linalg as la
import numpy as np

N = 25
L = 1.
dx = L/N
dx2 = dx**2
K1 = sparse.diags([1./dx2, -2/dx2, 1/dx2], [-1, 0, 1], shape=(N, N))
I = sparse.identity(N, format='dia')
K2 = sparse.kron(K1,I)+sparse.kron(I,K1)
I2 = sparse.kron(I,I)
K3 = sparse.kron(I,K2) + sparse.kron(K1,I2)

row  = np.array([N/2])
col  = np.array([0])
data = np.array([1/dx])
delta = sparse.coo_matrix((data, (row, col)), shape=(N,1)).tocsr()

source2 = sparse.kron(delta,delta)
source3 = sparse.kron(delta,source2)

phi = la.spsolve(K3, source3)
phi = phi.reshape((N,N,N))

x = np.linspace(0,L,N,endpoint=False)
y = np.linspace(0,L,N,endpoint=False)
z = np.linspace(0,L,N,endpoint=False)
xv, yv, zv = np.meshgrid(x,y,z)

from skimage import measure
verts, faces = measure.marching_cubes(phi, -.1)
#print verts
verts = verts * L / N # rescale verts
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts[faces])

#ax.set_xlim(0, 24)  # a = 6 (times two for 2nd ellipsoid)
#ax.set_ylim(0, 20)  # b = 10
#ax.set_zlim(0, 32)


Looks like a sphere. Cool.

Using scikit-image for finding the isosurface (constant potential surface). Ripped right from their examples.

Note that it needs to be rescaled, since verts is in the integer index format.

Seems to work. Takes a couple seconds to finish on my macbook pro at N=25.


NMR: a baby birds first fall from the nest

The idea is to attempt to use commodity hardware (standard hacker toolkit) to see an NMR signal

We bought 3 disk magnets from amazon. 1.5″x.25″.

Using our magnetometer made from a hall effect sensor, we’re getting that that field in the bore ~0.5T.

Going off a guess that disk magnets have fields approximately that of coils, a helmholtz coil arrangement was attempted with 3d printed parts. The fit of the magnets in there is really tight.

We saturated the sensor (which saturates at 0.3T), but with theoretically the thing is taking the dot product with it’s orientation with respect to the field. A perpendicular sensor does read approximately 0T. So we attenuated the field by tilting the sensor 45/2 degrees off of perpendicular and measured ~0.17T.

For the hydrogen nucleus, this is a resoannce freqeuncy of ~20Mhz.

We designed a helmholtz coil setup (coils that are 1 radius apart. They look unnaturally close to my eye for some reason, but that it just opinion.), Oversizing things that have to fit by 0.01inch seems enough for a snug fit. Your mileage may vary.

It is quite unclear to me the electrical characteristics these coils are going to have at the design frequency of ~20Mhz.

Here is a GUESS. I haven’t really done much hard work on this. But it’s probably right to order of magnitude. The B is fairly uniform in the helmholtz coil. Basically the same result for a solenoid of same size.

N is number of turns on one coil.

B = (4/5)^{3/2} \mu_0 N I / r

\Phi = L I

\Phi = 2 N A B

L =(4/5)^{3/2} \mu_0 N^2 2 \pi r

A = \pi r^2


Gives an inductance of 7uH, impedance of ~1kOhm at 20Mhz with 10 turns per coil.

This estimate was confirmed

Interesting Command Line Guys: Sox and Socat

Check out Socat. Let’s you pipe a file quickly over tcp and then grab it with netcat. Makes sense.

rtl_sdr -s 240000 -f 89500000 -g 20 - | tcc -lm -run minidemod-wfm.c \
    | sox -t raw -r 240000 -e unsigned -b 8 -c 1 - -t raw - rate 48000 \
    | mplayer -quiet -rawaudio samplesize=1:channels=1:rate=48000 -demuxer rawaudio -


Sox is a sound processing swiss army knife. Intriguing. Use it like image magick. Can quickly resample sound. That’s what the rate 48000 does. -r specifies that rtl-sdr is outputting 240000 samples per second. -b is that they are byte samples. -c 1 is one channel (as compared to stereo)


spectrogram is an interesting option. It can produce a spectrogram of the signal. Could use this as a waterfall plot.

synth is also interesting. Can be used to make waveforms.

sox −r 8000 −n output.wav synth 3 sine 300−3300

3 second sweep of a sine wave. 300Hz-3300Hz


A fun little man of rutherford scattering

So I coded up rutherford scattering in a real dumb way (you can significantly reduce your considerations by using symmetry and stuff).

I sort of monte carlo it with gaussian distributed initial conditions

import numpy as np
from scipy import integrate

z0 = -20.0
samples = 100
t = np.linspace(0,50,100)

def F(x):
    return .1 * x / np.linalg.norm(x)**1.5

def D(x,t):
    return np.concatenate((x[3:],F(x[0:3])))

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

def update(E):
    for i in range(samples):
        initxy = 5*np.random.randn(2)
        init = np.append(initxy,[z0, 0., 0., np.sqrt(E)])
        sol = integrate.odeint(D, init, t)
        ax.plot(sol[:,0], sol[:,1],sol[:,2])

The bundles that come off look pretty cool


Lots that one could do with this. Compare the outgoing distribution to the formula, Try to discern shape of other potentials. Do a little statistics to see if charge or whatever can be determined from the data.

Show center of mass scattering. Try 4 particle scattering.

I guess I’m trying to play around in some high energy concepts.


rtl_power & gnuradio

You can take a big honkin spectrum with rtl_power

rtl_power -f 118M:137M:8k -g 50 -i 10 -e 1h airband.csv

Frequency range and steps are self explanatory. Maximum steps is 2M ish (max sampling rate of SDR)

-i is integration time for each step (defualt is 10s)

-e is timeout time (will go forever if you don’t set)

-g is gain

Each line of the output is one window of sampling basically. The db values at the end are interpolating the frequencies in the window of the line.

Set -i and -e to the same number to just take one spectrum.  converts to a more sane format for my purposes, a list of frequencies and powers.

Then a short python script

import numpy as np
import sys
import matplotlib.pyplot as plt

filename = sys.argv[1]

freq, db = np.loadtxt(filename, delimiter=',', usecols=(0, 1), unpack=True)


that talks the csv filename as an argument will make a simple plot. You could also do this with plt.plotfile maybe? Something to investigate in the future. numpy.loadtxt is pretty clutch.


Or you can just spreadsheet it up

or use gnuplot (ick)




Also,  you can get a simple scope in gnuradio by just taking the rtl_sdr block and hook it right up to a scope block. Pretty Handy little guy.

Autocompleting gr_ ( double tab ) on my comp gives a list of interesting built ins:

gr_constellation_plot  gr_plot_psd_c          gr_spectrogram_plot_c

gr_filter_design       gr_plot_psd_f          gr_spectrogram_plot_f

gr_modtool             gr_plot_qt             gr_spectrogram_plot_i

gr_plot_char           gr_plot_short          gr_spectrogram_plot_s

gr_plot_const          gr_psd_plot_b          gr_time_plot_b

gr_plot_fft            gr_psd_plot_c          gr_time_plot_c

gr_plot_fft_c          gr_psd_plot_f          gr_time_plot_f

gr_plot_fft_f          gr_psd_plot_i          gr_time_plot_i

gr_plot_float          gr_psd_plot_s          gr_time_plot_s

gr_plot_int            gr_read_file_metadata  gr_time_raster_b

gr_plot_iq             gr_spectrogram_plot    gr_time_raster_f

gr_plot_psd            gr_spectrogram_plot_b 

Presumably these are for giving files to. The b, f, i, s must be the value type that was saved to file (byte char short float?)

For some reason I found gnuradio terrifying upon first inspection, but now it seems so friendly. What up wit dat?