imshow question

Dear All,

as i'm new to this list let me first introduce myself before posting a

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.

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,
                                        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
        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)
        rank = 1
    if (isinstance(input, (types.IntType, types.LongType, types.FloatType))):
        normalized = [input] * rank
        if isinstance(input, numarray.numarraycore.NumArray):
            normalized = list(input)
            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

    if PLOT_1D_GAUSS:

        # testing the Gaussiand Derivatives

        title('Gaussian Derivatives')

        scale = 9
        def _gDplot(n,ox,oy):
            imshow( gD(a,scale,(ox,oy)) )

        a = zeros( (64,64) ) * 1.0
        a[32,32] = 1.0
        title('Gaussian Derivatives')

        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
        scales = array([1,3,5,7,9,11,15,19, 25,31,41,51])
        plot( scales, map( time_gD, scales ), 'x' )
        xlabel('Scale (px)')
        ylabel('Time (s)')
        title('Gaussian Convolution Timings')