2011-11-16 129 views
1

我有一個圖像處理問題,我目前使用numpy和scipy在python中解決。簡而言之,我有一個想要應用許多局部收縮的圖像。我的原型代碼正在工作,最終的圖像看起來很棒。但是,處理時間已成爲我們應用程序中的嚴重瓶頸。你能幫我加快我的圖像處理代碼嗎?在python中加速插值

我試圖將我們的代碼歸結爲下面的'卡通'版本。剖析表明我正在花大部分時間進行插值。有沒有明顯的方法來加速執行?

import cProfile, pstats 
import numpy 
from scipy.ndimage import interpolation 

def get_centered_subimage(
    center_point, window_size, image): 
    x, y = numpy.round(center_point).astype(int) 
    xSl = slice(max(x-window_size-1, 0), x+window_size+2) 
    ySl = slice(max(y-window_size-1, 0), y+window_size+2) 
    subimage = image[xSl, ySl] 
    interpolation.shift(
     subimage, shift=(x, y)-center_point, output=subimage) 
    return subimage[1:-1, 1:-1] 

"""In real life, this is experimental data""" 
im = numpy.zeros((1000, 1000), dtype=float) 
"""In real life, this mask is a non-zero pattern""" 
window_radius = 10 
mask = numpy.zeros((2*window_radius+1, 2*window_radius+1), dtype=float) 
"""The x, y coordinates in the output image""" 
new_grid_x = numpy.linspace(0, im.shape[0]-1, 2*im.shape[0]) 
new_grid_y = numpy.linspace(0, im.shape[1]-1, 2*im.shape[1]) 


"""The grid we'll end up interpolating onto""" 
grid_step_x = new_grid_x[1] - new_grid_x[0] 
grid_step_y = new_grid_y[1] - new_grid_y[0] 
subgrid_radius = numpy.floor(
    (-1 + window_radius * 0.5/grid_step_x, 
    -1 + window_radius * 0.5/grid_step_y)) 
subgrid = (
    window_radius + 2 * grid_step_x * numpy.arange(
     -subgrid_radius[0], subgrid_radius[0] + 1), 
    window_radius + 2 * grid_step_y * numpy.arange(
     -subgrid_radius[1], subgrid_radius[1] + 1)) 
subgrid_points = ((2*subgrid_radius[0] + 1) * 
        (2*subgrid_radius[1] + 1)) 

"""The coordinates of the set of spots we we want to contract. In real 
life, this set is non-random:""" 
numpy.random.seed(0) 
num_points = 10000 
center_points = numpy.random.random(2*num_points).reshape(num_points, 2) 
center_points[:, 0] *= im.shape[0] 
center_points[:, 1] *= im.shape[1] 

"""The output image""" 
final_image = numpy.zeros(
    (new_grid_x.shape[0], new_grid_y.shape[0]), dtype=numpy.float) 

def profile_me(): 
    for m, cp in enumerate(center_points): 
     """Take an image centered on each illumination point""" 
     spot_image = get_centered_subimage(
      center_point=cp, window_size=window_radius, image=im) 
     if spot_image.shape != (2*window_radius+1, 2*window_radius+1): 
      continue #Skip to the next spot 
     """Mask the image""" 
     masked_image = mask * spot_image 
     """Resample the image""" 
     nearest_grid_index = numpy.round(
       (cp - (new_grid_x[0], new_grid_y[0]))/
       (grid_step_x, grid_step_y)) 
     nearest_grid_point = (
      (new_grid_x[0], new_grid_y[0]) + 
      (grid_step_x, grid_step_y) * nearest_grid_index) 
     new_coordinates = numpy.meshgrid(
      subgrid[0] + 2 * (nearest_grid_point[0] - cp[0]), 
      subgrid[1] + 2 * (nearest_grid_point[1] - cp[1])) 
     resampled_image = interpolation.map_coordinates(
      masked_image, 
      (new_coordinates[0].reshape(subgrid_points), 
      new_coordinates[1].reshape(subgrid_points)) 
      ).reshape(2*subgrid_radius[1]+1, 
         2*subgrid_radius[0]+1).T 
     """Add the recentered image back to the scan grid""" 
     final_image[ 
      nearest_grid_index[0]-subgrid_radius[0]: 
      nearest_grid_index[0]+subgrid_radius[0]+1, 
      nearest_grid_index[1]-subgrid_radius[1]: 
      nearest_grid_index[1]+subgrid_radius[1]+1, 
      ] += resampled_image 

cProfile.run('profile_me()', 'profile_results') 
p = pstats.Stats('profile_results') 
p.strip_dirs().sort_stats('cumulative').print_stats(10) 

的代碼做什麼含糊的解釋:

我們先從一個像素化的2D圖像,以及一組任意的(X,Y)在我們的形象加分不一般落在整數格。對於每個(x,y)點,我想通過一個小屏蔽中心精確地在該點上乘以圖像。接下來,我們將有限量的遮罩區域縮小/展開,最後將這個處理後的子圖像添加到可能與原始圖像不具有相同像素大小的最終圖像。 (不是我最好的解釋,好吧)。

+1

作爲第一次切割,你可能想* [給這個嘗試](http://stackoverflow.com/questions/4295799/how-to-improve-performance-of-this-code/4299378#4299378) *。 –

回答

3

我敢肯定,正如你所說的,大部分計算時間發生在interpolate.map_coordinates(…)之間,在center_points的每次迭代中調用一次,這裏是10,000次。一般來說,使用numpy/scipy堆棧時,您希望重複執行大型數組的任務在本地Numpy/Scipy函數中發生 - 即在同類數據的C循環中 - 而不是在Python中顯式執行。

一個策略可能加快插值,但也將增加使用的內存量,方法是:

  • 首先,獲取在三維陣列中的所有子圖像(這裏命名masked_image)( window_radius x window_radius x center_points.size
  • 編寫一個ufunc(讀取它很有用),它使用numpy.frompyfunc包裝每個子圖像必須完成的工作,該工作應該返回另一個三維數組(subgrid_radius[0] x subgrid_radius[1] x center_points.size)。簡而言之,這會創建一個向量化的python函數版本,它可以在數組上按元素進行廣播。
  • 通過在第三維上求和來構建最終圖像。

希望讓你更接近你的目標!