1 import numpy as np
   2 cimport numpy as np
   3 
   4 cdef void _advance(float vx, float vy,
   5         int* x, int* y, float*fx, float*fy, int w, int h):
   6     cdef float tx, ty
   7     if vx>=0:
   8         tx = (1-fx[0])/vx
   9     else:
  10         tx = -fx[0]/vx
  11     if vy>=0:
  12         ty = (1-fy[0])/vy
  13     else:
  14         ty = -fy[0]/vy
  15     if tx<ty:
  16         if vx>=0:
  17             x[0]+=1
  18             fx[0]=0
  19         else:
  20             x[0]-=1
  21             fx[0]=1
  22         fy[0]+=tx*vy
  23     else:
  24         if vy>=0:
  25             y[0]+=1
  26             fy[0]=0
  27         else:
  28             y[0]-=1
  29             fy[0]=1
  30         fx[0]+=ty*vx
  31     if x[0]>=w:
  32         x[0]=w-1 # FIXME: other boundary conditions?
  33     if x[0]<0:
  34         x[0]=0 # FIXME: other boundary conditions?
  35     if y[0]<0:
  36         y[0]=0 # FIXME: other boundary conditions?
  37     if y[0]>=h:
  38         y[0]=h-1 # FIXME: other boundary conditions?
  39 
  40 
  41 #np.ndarray[float, ndim=2] 
  42 def line_integral_convolution(
  43         np.ndarray[float, ndim=3] vectors,
  44         np.ndarray[float, ndim=2] texture,
  45         np.ndarray[float, ndim=1] kernel):
  46     cdef int i,j,k,x,y
  47     cdef int h,w,kernellen
  48     cdef int t
  49     cdef float fx, fy, tx, ty
  50     cdef np.ndarray[float, ndim=2] result
  51 
  52     h = vectors.shape[0]
  53     w = vectors.shape[1]
  54     t = vectors.shape[2]
  55     kernellen = kernel.shape[0]
  56     if t!=2:
  57         raise ValueError("Vectors must have two components (not %d)" % t)
  58     result = np.zeros((h,w),dtype=np.float32)
  59 
  60     for i in range(h):
  61         for j in range(w):
  62             x = j
  63             y = i
  64             fx = 0.5
  65             fy = 0.5
  66             
  67             k = kernellen//2
  68             #print i, j, k, x, y
  69             result[i,j] += kernel[k]*texture[x,y]
  70             while k<kernellen-1:
  71                 _advance(vectors[y,x,0],vectors[y,x,1],
  72                         &x, &y, &fx, &fy, w, h)
  73                 k+=1
  74                 #print i, j, k, x, y
  75                 result[i,j] += kernel[k]*texture[x,y]
  76 
  77             x = j
  78             y = i
  79             fx = 0.5
  80             fy = 0.5
  81             
  82             while k>0:
  83                 _advance(-vectors[y,x,0],-vectors[y,x,1],
  84                         &x, &y, &fx, &fy, w, h)
  85                 k-=1
  86                 #print i, j, k, x, y
  87                 result[i,j] += kernel[k]*texture[x,y]
  88                     
  89     return result