Dear All,
as i'm new to this list let me first introduce myself before posting a
question.
I'm working in the field of image processing and computer vision for a long
time now and i am a long time matlab user. Up to now i didn't like any of the
matlab clones but it seems that numarray+matplotlib might be a winner here.
I started using numarray+matplotlib by coding the Gaussian derivative
convolutions. Calculating a derivative (by convolving with the derivative of
the Gaussian function) will lead to an image(array) with both positive and
negative values. Although imshow should deal with that (as far as i
understand the code: beware i am a Python beginner) the display turns black
for all derivatives (except for the 'zero order' derivative of course).
What is going on? Is it me (probably) or is it imshow?
Attached you find the code.
Regards
Rein van den Boomgaard
University of Amsterdam
The Netherlands
# ==========================================================
# Gaussian Derivatives
···
#
# Copyright (c) Rein van den Boomgaard
# Vision Consultancy & Training
# rein@...237...
#
# ==========================================================
import types
from numarray import *
from numarray.nd_image import *
def gD( image, scales, orders, nscales=4, mode='reflect' ):
scales = _normalize_sequence(scales,image)
orders = _normalize_sequence(orders,image)
result = image.copy()
ndims = len( image.shape )
for i in range(0,ndims):
N = ceil(scales[i]*nscales)
x = arange( -N, N )
w = GaussianDerivativeFunction( x,
scales[i],
orders[i] )
result = convolve1d( result, w, i, mode = mode )
return( result )
def HermiteH( n, x ):
"""
Hermite Poynomial
"""
if n==0:
return 0*x + 1
elif n==1:
return 2*x
else:
return 2*x*HermiteH( n-1, x ) - 2*(n-1)*HermiteH(n-2,x)
def GaussianDerivativeFunction( x, scale, order ):
scale = float(scale)
g = 1.0 / ( scale * sqrt(2*pi) ) * exp( - x**2 / (2*scale**2) )
return (-1.0 / (scale*sqrt(2.0)))**order * \
HermiteH( order, x / (scale * sqrt(2)) ) * g
# the _normalize_sequence function is copied from numarray
# because i cant import it...???
def _normalize_sequence(input, array):
"""If input is a scalar, create a sequence of length equal to the
rank of array by duplicating the input. If input is a sequence,
check if its length is equalt to the lenght of array.
"""
if isinstance(array, ArrayType):
rank = len(array.shape)
else:
rank = 1
if (isinstance(input, (types.IntType, types.LongType, types.FloatType))):
normalized = [input] * rank
else:
if isinstance(input, numarray.numarraycore.NumArray):
normalized = list(input)
else:
normalized = input
if len(normalized) != rank:
error = "sequence argument must have length equal to input rank"
raise RuntimeError, error
return normalized
if __name__ == "__main__":
from matplotlib import use, interactive
from matplotlib.matlab import *
PLOT_1D_GAUSS = 1
PLOT_GAUSSIAN_DERIVATIVES = 1
PLOT_TIMING = 1
if PLOT_1D_GAUSS:
# testing the Gaussiand Derivatives
x=arange(-5,5,.01)
figure(1)
clf()
plot(x,GaussianDerivativeFunction(x,1,0))
plot(x,GaussianDerivativeFunction(x,1,1))
plot(x,GaussianDerivativeFunction(x,1,2))
plot(x,GaussianDerivativeFunction(x,1,3))
plot(x,GaussianDerivativeFunction(x,1,4))
xlabel(r'$x$')
ylabel(r'$G_n(x)$')
title('Gaussian Derivatives')
if PLOT_GAUSSIAN_DERIVATIVES:
scale = 9
def _gDplot(n,ox,oy):
subplot(4,4,n)
axis('off')
imshow( gD(a,scale,(ox,oy)) )
a = zeros( (64,64) ) * 1.0
a[32,32] = 1.0
figure(2)
axis('off')
title('Gaussian Derivatives')
_gDplot(1,0,0)
title(r'$G$')
_gDplot(5,1,0)
title(r'$G_x$')
_gDplot(6,0,1)
title(r'$G_y$')
_gDplot(9,2,0)
title(r'$G_{xx}$')
_gDplot(10,1,1)
title(r'$G_{xy}$')
_gDplot(11,0,2)
title(r'$G_{yy}$')
_gDplot(13,3,0)
title(r'$G_{xxx}$')
_gDplot(14,2,1)
title(r'$G_{xxy}$')
_gDplot(15,1,2)
title(r'$G_{xyy}$')
_gDplot(16,0,3)
title(r'$G_{yyy}$')
if PLOT_TIMING:
from timeit import Timer
niter = 5
def time_gD(scale):
t = Timer( "b=gD(a,%d, 0)" % scale, "from __main__ import gD,a" )
timing = t.timeit(number=niter )
print( timing/niter )
return( timing/niter )
a = zeros( (256,256) ) * 1.0
a[128,128] = 1.0
figure(3)
scales = array([1,3,5,7,9,11,15,19, 25,31,41,51])
print(type(scales))
plot( scales, map( time_gD, scales ), 'x' )
xlabel('Scale (px)')
ylabel('Time (s)')
title('Gaussian Convolution Timings')