CODEWITHSUNDEEP
pythonnumpyscipyleast-squares
Craig
Craig

Reputation: 41

Scipy leastsq() function overhead

I am working on an image analysis program and I have narrowed down my bottleneck to attempts to fit a 2D gaussian to a small window (20x20) pixels many times. 90% of the execution time is spent in this code.

I am using the code given in the scipy cookbook for this problem:

   def gaussian(height, center_x, center_y, width_x, width_y):
           """Returns a gaussian function with the given parameters"""
        width_x = float(width_x)
        width_y = float(width_y)
        return lambda x,y: height*exp(
                    -(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)


def moments(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution by calculating its
    moments """
    total = data.sum()
    X, Y = indices(data.shape)
    x = (X*data).sum()/total
    y = (Y*data).sum()/total
    col = data[:, int(y)]
    width_x = sqrt(abs((arange(col.size)-y)**2*col).sum()/col.sum())
    row = data[int(x), :]
    width_y = sqrt(abs((arange(row.size)-x)**2*row).sum()/row.sum())
    height = data.max()
    return height, x, y, width_x, width_y


def fitgaussian(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution found by a fit"""
    params = moments(data)
    errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) -
                                 data)
    p, success = optimize.leastsq(errorfunction, params, maxfev=50, ftol=1.49012e-05)
    return p

I was able to cut the execution time in half by combining the errorfunction() and gaussian() functions so every time leastsq() calls errorfunction() there is one function call instead of two.

This leads me to believe that most of the remaining execution time is spent tied up in function call overhead as the leastsq() algorithm calls errorfunction().

Is there any way to reduce this function call overhead? I am at a loss since leastsq() takes a function as an input.

I apologize in advance if my description is confusing, I am a mechanical engineer by training and I'm learning Python as I go. Please let me know if there is any other information that would be helpful.

Upvotes: 4

Views: 569

Answers (1)

mtadd
mtadd

Reputation: 2555

Since exp is monotonic, you could use the logarithm of the gaussian as your error function, e.g.

def log_gaussian(height, center_x, center_y, width_x, width_y):
    """Returns a gaussian function with the given parameters"""
    width_x = float(width_x)
    width_y = float(width_y)
    log_height = log(height)
    return lambda x,y: (log_height - 
             (((center_x-x)/width_x)**2 - ((center_y-y)/width_y)**2)/2)

This will result in 1 call to log per iteration, rather than 1 call to exp per dataset row each iteration.

Upvotes: 1

Related Questions