Webcam image processing

This script acquires a digital camera image from a Logitech webcam and processes it. I have been using it for online monitoring of an atomic hydrogen beam on a beam-dump quartz plate. To get repeatable images, it is helpful to set the webcam's gain and exposure.

	# Display webcam image from beamline monitor
	# version 3, 2014-04-25
	# Amar
	import numpy as np
	import cv2
	import matplotlib.animation as animation
	import matplotlib.pyplot as plt
	import matplotlib.cm as cm
	from matplotlib.figure import Figure
	import matplotlib.gridspec as gridspec
	class Camera:
	def __init__(self,channel=0):
	self.capture = cv2.VideoCapture(channel)
	self.success,self.image = self.capture.read()
	print "Beginning acquisition ...",
	def acquire(self):
	# Acquire image
	self.success,self.image = self.capture.read()
	print ".",
	def close(self):
	if self.capture.isOpened(): self.capture.release()
	print "released camera"
	camera = Camera()
	camera.capture.set(15,-3)
	camera.capture.set(16,0.01)
	"""
	Camera parameters:
	--------------------
	1 CV_CAP_PROP_POS_MSEC Current position of the video file in milliseconds.
	2 CV_CAP_PROP_POS_FRAMES 0-based index of the frame to be decoded/captured next.
	3 CV_CAP_PROP_POS_AVI_RATIO Relative position of the video file
	4 CV_CAP_PROP_FRAME_WIDTH Width of the frames in the video stream.
	5 CV_CAP_PROP_FRAME_HEIGHT Height of the frames in the video stream.
	6 CV_CAP_PROP_FPS Frame rate.
	7 CV_CAP_PROP_FOURCC 4-character code of codec.
	8 CV_CAP_PROP_FRAME_COUNT Number of frames in the video file.
	9 CV_CAP_PROP_FORMAT Format of the Mat objects returned by retrieve() .
	10 CV_CAP_PROP_MODE Backend-specific value indicating the current capture mode.
	11 CV_CAP_PROP_BRIGHTNESS Brightness of the image (only for cameras).
	12 CV_CAP_PROP_CONTRAST Contrast of the image (only for cameras).
	13 CV_CAP_PROP_SATURATION Saturation of the image (only for cameras).
	14 CV_CAP_PROP_HUE Hue of the image (only for cameras).
	15 CV_CAP_PROP_GAIN Gain of the image (only for cameras).
	16 CV_CAP_PROP_EXPOSURE Exposure (only for cameras).
	17 CV_CAP_PROP_CONVERT_RGB Boolean flags indicating whether images should be converted to RGB.
	18 CV_CAP_PROP_WHITE_BALANCE Currently unsupported
	19 CV_CAP_PROP_RECTIFICATION Rectification flag for stereo cameras (note: only supported by DC1394 v 2.x backend currently)
	"""
	# Setting up the plot
	fig = plt.figure()
	aspectRatio = 6
	fig.canvas.set_window_title('beam camera')
	gs = gridspec.GridSpec(aspectRatio,aspectRatio)
	ax1 = plt.subplot(gs[0,:-1])
	ax2 = plt.subplot(gs[1:,:-1])
	ax3 = plt.subplot(gs[1:,-1])
	plt.subplots_adjust(left=0.10,bottom=0.10,wspace=0.02, hspace=0.02)
	ax1.xaxis.set_visible(False)
	ax1.yaxis.set_visible(False)
	ax3.xaxis.set_visible(False)
	ax3.yaxis.set_visible(False)
	ax1.set_xlim(0,640)
	ax3.set_ylim(480,0)
	### Initial plot ###
	image = 0.5*camera.image[:,:,0] + 0.5*camera.image[:,:,1]
	oldImage = camera.image
	i,j = np.unravel_index(image.argmax(), image.shape)
	print i,j
	x = image[i] #x slice along maximum
	y = image.T[j] #y slice along maximum
	xslice, = ax1.plot(range(640),x,'b',linewidth=3,alpha=0.7)
	im = ax2.imshow(image)
	yslice, = ax3.plot(y[::-1],range(480,0,-1),'b',linewidth=3,alpha=0.7)
	text_template = 'x0 = %i, y0 = %i'
	text = ax2.text(300,40,text_template%(j,i),family='monospace',fontsize=16,weight='heavy',color='white')
	#### Animation routine ###
	def updatefig(*args):
	camera.acquire()
	image = 0.5*camera.image[:,:,0] + 0.5*camera.image[:,:,1]
	if camera.success:
	image = 0.5*camera.image[:,:,0] + 0.5*camera.image[:,:,1]
	oldImage = camera.image
	i,j = np.unravel_index(image.argmax(), image.shape)
	print i,j
	x = image[i] #x slice along maximum
	y = image.T[j] #y slice along maximum
	im.set_array(image)
	xslice.set_ydata(x)
	ax1.set_ylim(0,x.max())
	yslice.set_xdata(y[::-1])
	ax3.set_xlim(0,y.max())
	text.set_text(text_template%(j,i))
	else: print "error"
	return (im,) + (xslice,) + (yslice,) + (text,)
	def init(): #Init only required for blitting to give a clean slate.
	im.set_array(np.zeros(oldImage.shape))
	xslice.set_ydata(np.zeros((640,)))
	yslice.set_xdata(np.zeros((480,)))
	text.set_text("")
	return (im,) + (xslice,) + (yslice,) + (text,)
	ani = animation.FuncAnimation(fig, updatefig, init_func=init, interval=80, blit=True)
	plt.ioff()
	plt.show()
	camera.close()

view raw beamCamera.py hosted with ❤ by GitHub

Oscillator animation

An animated simple harmonic oscillator, integrated using the Verlet algorithm. You can see that the Verlet algorithm conserves energy, but only in an average sense.

	# Animation of harmonic oscillator, with Verlet integration
	from __future__ import division
	import numpy as np
	from numpy import sin, cos
	from scipy.constants import c,e,h,hbar,u,pi
	import matplotlib.animation as animation
	delta = 1e-2 # spatial discretization
	x = np.arange(-2,2,delta)
	fig, ax = plt.subplots()
	# ***** ANIMATION ***********
	def V(x): return 0.5*x**2 # m=1, omega=1
	def a(x,d=delta): return -( V(x+d) - V(x) )/d
	def energy(x,v): return 0.5*v**2 + V(x)
	line, = ax.plot(x,V(x),lw=2)
	x0,v0 = np.random.uniform(-1,1,size=2)
	blob, = ax.plot(x0,V(x0),'ro',ms=16)
	xi,vi = x0,v0
	text_template = 'energy = %.2g'
	text = ax.text(0.25,0.5,text_template%(energy(x0,v0)), transform=ax.transAxes, family='monospace',fontsize=20,weight='heavy')
	def animate(i,h=1e-1):
	global xi, vi, a
	t = i*h
	## Verlet algorithm
	xold = xi
	xi += vi*h + 0.5*a(xi)*h**2
	vi += 0.5*( a(xold) + a(xi) )*h
	blob.set_xdata( xi )
	blob.set_ydata( V(xi) ) # update the data
	text.set_text( text_template%(energy(xi,vi)) )
	return (blob,) + (text,)
	def init(): #Init only required for blitting to give a clean slate.
	blob.set_ydata(np.ma.array(x, mask=True))
	text.set_text(" ")
	return (blob,) + (text,)
	ani = animation.FuncAnimation(fig, animate, range(1000), init_func=init,
	interval=0.2, blit=True)
	plt.show()

view raw oscillatorAnimation.py hosted with ❤ by GitHub

Bringing fourier transforms to life

This script animates the reconstruction of a random ellipse from its fourier coefficients. It is a lot of fun to watch the reconstruction, especially for large aspect ratio ellipses.

	# Fourier reconstruction of random ellipse
	# Amar Vutha
	# 2014-04-03
	import pylab as plt
	import numpy as np
	from numpy import pi
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	# Generate random ellipse
	t = np.arange(0,2*pi,1e-4)
	r = np.zeros(len(t))
	a,b = np.random.random(), np.random.random()
	theta = t.copy()
	r = a*b/np.sqrt( ( a*np.sin(theta) )**2 + ( b*np.cos(theta) )**2 )
	fig = plt.figure()
	ax = plt.subplot(111,polar=True)
	ax.set_rlim(0,1)
	ax.plot(theta,r,'r--',lw=3,alpha=0.5)
	line, = ax.plot(theta,r,lw=2)
	linebg, = ax.plot(theta,np.ones(len(theta)),'g--',lw=3,alpha=0.2)
	### Fourier transform #####
	phi = t.copy()
	Nmax = 20
	A,B = [],[]
	dp = phi[1] - phi[0]
	A_DC = np.sum(r) * dp/(2*pi)
	A,B = np.zeros( (Nmax) ), np.zeros( (Nmax) )
	for n in range(1,Nmax):
	A[n] = np.sum(r*np.cos(n*phi)) * dp/pi
	B[n] = np.sum(r*np.sin(n*phi)) * dp/pi
	### Reconstruction ######
	def reconstruct(M):
	R = A_DC
	for m in range(1,M):
	R += A[m]* np.cos(m*phi) + B[m]* np.sin(m*phi)
	return R
	def animate(i):
	line.set_ydata(reconstruct(i)) # update the data
	linebg.set_ydata(np.cos(i*theta)) # shows the shape corresponding to the fourier coefficient
	return (line,) + (linebg,)
	#Init only required for blitting to give a clean slate.
	def init():
	line.set_ydata(np.ma.array(phi, mask=True))
	linebg.set_ydata(np.ma.array(phi, mask=True))
	return (line,) + (linebg,)
	ani = animation.FuncAnimation(fig, animate, np.arange(1, Nmax), init_func=init,
	interval=500, blit=True)
	plt.show()

view raw fourier-ellipse.py hosted with ❤ by GitHub