[Matplotlib-users] imshow question

[Rein v. d. B. (UvA) <re...@sc...>]

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
#               re...@vi...
# 
# ==========================================================
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')