Function iapconv

  • 2D or 3D Periodic Convolution.

Synopse

1D,2D or 3D Periodic convolution.

  • g = iapconv(f, h)
    • g: Image.
    • f: Image. Input image.
    • h: Image. PSF (point spread function), or kernel. The origin is at the array center.

Description

  • Perform a 1D, 2D or 3D discrete periodic convolution. The kernel origin is at the center of image h. Both image and kernel are periodic with same period. Usually the kernel h is smaller than the image f, so h is padded with zero until the size of f.
01. def iapconv4(f,h):
02.    import numpy as np
03. 
04.    h_ind=np.nonzero(h)
05.    f_ind=np.nonzero(f)
06.    if len(h_ind[0])>len(f_ind[0]):
07.        h,    f    = f,    h
08.        h_ind,f_ind= f_ind,h_ind
09. 
10.    gs = np.maximum(np.array(f.shape),np.array(h.shape))
11.    if (f.dtype == 'complex') or (h.dtype == 'complex'):
12.        g = np.zeros(gs,dtype='complex')
13.    else:
14.        g = np.zeros(gs)
15. 
16.    f1 = g.copy()
17.    f1[f_ind]=f[f_ind]
18. 
19.    if f.ndim == 1:
20.      (W,) = gs
21.      col = np.arange(W)
22.      for cc in h_ind[0]:
23.         g[:] += f1[(col-cc)%W] * h[cc]
24. 
25.    elif f.ndim == 2:
26.      H,W = gs
27.      row,col = np.indices(gs)
28.      for rr,cc in np.transpose(h_ind):
29.         g[:] += f1[(row-rr)%H, (col-cc)%W] * h[rr,cc]
30. 
31.    else:
32.      Z,H,W = gs
33.      d,row,col = np.indices(gs)
34.      for dd,rr,cc in np.transpose(h_ind):
35.         g[:] += f1[(d-dd)%Z, (row-rr)%H, (col-cc)%W] * h[dd,rr,cc]
36.    return g
01. from numpy import *
02. 
03. def iapconv(f, h):
04. 
05.     f, h = asarray(f), asarray(h,dtype=float)
06.     #faux, haux = ravel(f), ravel(h)
07. 
08.     s = f.shape
09.     # Checking the shape of the image F.
10.     if len(f.shape) == 1:
11.         f = f[newaxis,newaxis,:]
12.     elif len(f.shape) == 2:
13.         f = f[newaxis,:,:]
14. 
15.     # Checking the shape of the image H.
16.     if len(h.shape) == 1:
17.         h = h[newaxis,newaxis,:]
18.     elif len(h.shape) == 2:
19.         h = h[newaxis,:,:]
20. 
21.     # Getting the dimensions of images F and H
22.     (fslices, frows, fcols) = f.shape
23.     (hslices, hrows, hcols) = h.shape
24. 
25.     ds1 = int((hslices-1)/2.)
26.     ds2 = hslices - ds1
27. 
28.     dr1 = int((hrows-1)/2.)
29.     dr2 = hrows - dr1
30. 
31.     dc1 = int((hcols-1)/2.)
32.     dc2 = hcols - dc1
33. 
34.     p = concatenate((concatenate((f[-ds2+1::,:,:], f)), f[0:ds1,:,:]))
35.     p = concatenate((concatenate((p[:,-dr2+1::,:], p), 1), p[:,0:dr1,:]), 1)
36.     p = concatenate((concatenate((p[:,:,-dc2+1::], p), 2), p[:,:,0:dc1]), 2)
37. 
38.     g = zeros((fslices,frows,fcols))
39. 
40.     for i in range(hslices):
41.         for j in range(hrows):
42.             for k in range(hcols):
43.                 hw = h[hslices-i-1, hrows-j-1, hcols-k-1]
44.                 if (hw):
45.                     g = g + h[hslices-i-1,hrows-j-1,hcols-k-1] * p[i:fslices + i,j:frows + j, k:fcols + k]
46.     g.shape = s
47.     return g

Examples

Numerical Example 1D

01. import ia636 as ia
02. 
03. f = array([1,0,0,0,0,1,0,0,0])
04. print "\n Image (f):"
05. print f
06. 
07. h = array([1,2,3])
08. print "\n Image Kernel (h):"
09. print h
10. 
11. result = ia.iapconv(f,h)
12. print "\n Image Output - (g):"
13. print result
Image (f):
[1 0 0 0 0 1 0 0 0]

 Image Kernel (h):
[1 2 3]

 Image Output - (g):
[ 2.  3.  0.  0.  1.  2.  3.  0.  1.]

Numerical Example 2D

01. f = array([[1,0,0,0,0,0,0,0,0],
02.            [0,0,0,0,0,0,0,0,0],
03.            [0,0,0,1,0,0,0,0,0],
04.            [0,0,0,0,0,0,0,0,1],
05.            [0,0,0,0,0,0,0,0,0]])
06. print "\n Image (f):"
07. print f
08. 
09. h = array([[1,2,3],
10.            [4,5,6]])
11. print "\n Image Kernel (h):"
12. print h
13. 
14. result = ia.iapconv(f,h)
15. print "\n Image Output - (G):"
16. print result
Image (f):
[[1 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 1 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 0 0]]

 Image Kernel (h):
[[1 2 3]
 [4 5 6]]

 Image Output - (G):
[[ 2.  3.  0.  0.  0.  0.  0.  0.  1.]
 [ 5.  6.  0.  0.  0.  0.  0.  0.  4.]
 [ 0.  0.  1.  2.  3.  0.  0.  0.  0.]
 [ 3.  0.  4.  5.  6.  0.  0.  1.  2.]
 [ 6.  0.  0.  0.  0.  0.  0.  4.  5.]]

Numerical Example 3D

01. f = zeros((3,3,3))
02. f[0,1,1] = 1
03. f[1,1,1] = 1
04. f[2,1,1] = 1
05. 
06. print "\n Image Original (F): "
07. print f
08. 
09. h = array([[[ 1,  2, 3 ],[ 3,  4, 5 ], [ 5,  6, 7 ]],
10.           [ [ 8,  9, 10],[11, 12, 13], [14, 15, 16]],
11.           [ [17, 18, 19],[20, 21, 22], [23, 24, 25]]])
12. 
13. print "\n Image Kernel (H): "
14. print h
15. 
16. result = ia.iapconv(f,h)
17. print "\n Image Output - (G): "
18. print result
Image Original (F): 
[[[ 0.  0.  0.]
  [ 0.  1.  0.]
  [ 0.  0.  0.]]

 [[ 0.  0.  0.]
  [ 0.  1.  0.]
  [ 0.  0.  0.]]

 [[ 0.  0.  0.]
  [ 0.  1.  0.]
  [ 0.  0.  0.]]]

 Image Kernel (H): 
[[[ 1  2  3]
  [ 3  4  5]
  [ 5  6  7]]

 [[ 8  9 10]
  [11 12 13]
  [14 15 16]]

 [[17 18 19]
  [20 21 22]
  [23 24 25]]]

 Image Output - (G): 
[[[ 26.  29.  32.]
  [ 34.  37.  40.]
  [ 42.  45.  48.]]

 [[ 26.  29.  32.]
  [ 34.  37.  40.]
  [ 42.  45.  48.]]

 [[ 26.  29.  32.]
  [ 34.  37.  40.]
  [ 42.  45.  48.]]]

Example with Image 2D

01. f = adread('cameraman.pgm')
02. adshow(f, title = 'a) - Original Image')
03. h = [[-1,-1,-1],[0,0,0],[1,1,1]]
04. g = ia.iapconv(f,h)
05. print "\nPrewitt´s Mask"
06. print h
07. 
08. gn = ia.ianormalize(g, [0,255])
09. adshow(gn, title = 'b) Prewitt´s Mask filtering')
10. 
11. adshow(ia.ianormalize(abs(g)), title = 'c) absolute of Prewitt´s Mask filtering')
Prewitt´s Mask
[[-1, -1, -1], [0, 0, 0], [1, 1, 1]]

a) - Original Image

b) Prewitt´s Mask filtering

c) absolute of Prewitt´s Mask filtering

Example with Image 3D

01. import dicom
02. import numpy as np
03. 
04. filename = 'PHILIPS/DICOM/IM_0007'
05. dataset = dicom.read_file(find_attachment_file(filename))
06. 
07. linhas = dataset.Rows
08. colunas = dataset.Columns
09. fatias = dataset.NumberofFrames
10. 
11. data = dataset.pixel_array.astype(float32)
12. data = data[90:96,:,:]
13. 
14. h = array([[[-1, -1, -1 ],
15.             [ 0,  0,  0],
16.             [ 1,  1,  1 ]],
17.           [ [-1, -1, -1 ],
18.             [ 0,  0,  0],
19.             [ 1,  1,  1]],
20.           [ [-1, -1, -1 ],
21.             [ 0,  0,  0],
22.             [ 1,  1,  1 ]]])
23. 
24. g = ia.iapconv(data,h)
25. 
26. gdata = ia.ianormalize(data).astype(uint8)
27. gn = ia.ianormalize(g).astype(uint8)
28. gn1 = ia.ianormalize(abs(g)).astype(uint8)
29. 
30. adshow(ia.iamosaic(gdata,3), title ='Figure a) - Original Image')
31. adshow(ia.iamosaic(gn,   3), title ='Figure b) - 3D Periodic Convolution. - Prewitt´s Mask')
32. adshow(ia.iamosaic(gn1,  3), title ='Figure c) - 3D Periodic Convolution. - Prewitt´s Mask (Negative)')

Figure a) - Original Image

Figure b) - 3D Periodic Convolution. - Prewitt´s Mask

Figure c) - 3D Periodic Convolution. - Prewitt´s Mask (Negative)

Equation

Limitations

Both image and kernel are internally converted to double.

See Also

Contributions

  • Danilo Rodrigues Pereira, 1st semester 2011