Thursday, August 25, 2011

How to use ginput

ginput is a function that enables you to select points from a figure using the mouse. This post is a simple example on how to use it.
from pylab import plot, ginput, show, axis

axis([-1, 1, -1, 1])
print "Please click three times"
pts = ginput(3) # it will wait for three clicks
print "The point selected are"
print pts # ginput returns points as tuples
x=map(lambda x: x[0],pts) # map applies the function passed as 
y=map(lambda x: x[1],pts) # first parameter to each element of pts
axis([-1, 1, -1, 1])
And after three clicks this is the result on the console:
Please click three times
The point selected are
[(-0.77468982630272942, -0.3418367346938776), (0.11464019851116632, -0.21428571428571436), (0.420347394540943, 0.55612244897959173)]
and this is the figure generated:

Wednesday, August 10, 2011

Applying a Moebius transformation to an image

The Moebius transformation is defined as
where z is a complex variable different than -d/c. And a, b, c, d are complex numbers.
In this post we will see how to apply the Moebius transform to an image.
Given a set of three distinct points on the complex plane z1, z2, z3 and a second set of distinct points w1, w2, w3, there exists precisely one Moebius transformation f(z) which maps the zs to the ws. So, in the first step we have to determine f(z) from the given sets of points. We will use the explicit determinant formula to compute the coefficients:
from pylab import *
from numpy import *

zp=[157+148j, 78+149j, 54+143j]; # (zs) the complex point zp[i]
wa=[147+143j, 78+140j, 54+143j]; # (ws) will be in wa[i]

# transformation parameters
a = linalg.det([[zp[0]*wa[0], wa[0], 1], 
                [zp[1]*wa[1], wa[1], 1], 
                [zp[2]*wa[2], wa[2], 1]]);

b = linalg.det([[zp[0]*wa[0], wa[0], wa[0]], 
                [zp[1]*wa[1], wa[1], wa[1]], 
                [zp[2]*wa[2], wa[2], wa[2]]]);         

c = linalg.det([[zp[0], wa[0], 1], 
                [zp[1], wa[1], 1], 
                [zp[2], wa[2], 1]]);

d = linalg.det([[zp[0]*wa[0], zp[0], 1], 
                [zp[1]*wa[1], zp[1], 1], 
                [zp[2]*wa[2], zp[2], 1]]);
Now we can apply the transformation to the pixel coordinates.
img = fliplr(imread('mondrian.jpg')) # load an image

r = ones((500,500,3),dtype=uint8)*255 # empty-white image
for i in range(img.shape[0]):
 for j in range(img.shape[1]):
  z = complex(i,j)
  # transformation applied to the pixels coordinates
  w = (a*z+b)/(c*z+d)
  r[int(real(w)),int(imag(w)),:] = img[i,j,:] # copy of the pixel

title('Original Mondrian')
title('Mondrian after the transformation')
And this is the result.
This is another image obtained changing the vectors zp and wa.
As we can see, the result of the transformation has some empty areas that we should fill using interpolation. So let's look to another way to implement a geometric transformation on an image. The following code uses the scipy.ndimage.geometric_transform to implement the inverse Moebius transform:
from scipy.ndimage import geometric_transform

def shift_func(coords):
 """ Define the moebius transformation, though backwards """
 #turn the first two coordinates into an imaginary number
 z = coords[0] + 1j*coords[1]
 w = (d*z-b)/(-c*z+a) #the inverse mobius transform
 #take the color along for the ride
 return real(w),imag(w),coords[2]
r = geometric_transform(img,shift_func,cval=255,output_shape=(450,350,3)) 
It gives the following result:
We can see that the geometric_transform provides the pixel interpolation automatically.

Wednesday, August 3, 2011

How to plot the frequency spectrum with scipy

Spectrum analysis is the process of determining the frequency domain representation of a time domain signal and most commonly employs the Fourier transform. The Discrete Fourier Transform (DFT) is used to determine the frequency content of signals and the Fast Fourier Transform (FFT) is an efficient method for calculating the DFT. Scipy implements FFT and in this post we will see a simple example of spectrum analysis:
from numpy import sin, linspace, pi
from pylab import plot, show, title, xlabel, ylabel, subplot
from scipy import fft, arange

def plotSpectrum(y,Fs):
 Plots a Single-Sided Amplitude Spectrum of y(t)
 n = len(y) # length of the signal
 k = arange(n)
 T = n/Fs
 frq = k/T # two sides frequency range
 frq = frq[range(n/2)] # one side frequency range

 Y = fft(y)/n # fft computing and normalization
 Y = Y[range(n/2)]
 plot(frq,abs(Y),'r') # plotting the spectrum
 xlabel('Freq (Hz)')

Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = arange(0,1,Ts) # time vector

ff = 5;   # frequency of the signal
y = sin(2*pi*ff*t)

The program shows the following figure, on top we have a plot of the signal and on the bottom the frequency spectrum.