iapconv2 - Periodic convolution, kernel origin at array origin

Synopse

1D, 2D or 3D Periodic convolution. (kernel origin at array origin)

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

Description

Perform a 1D, 2D or 3D discrete periodic convolution. The kernel origin is at the origin 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. Supports complex images.

01. def iapconv2(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

Examples

Numerical Example 1D

01. import ia636 as ia
02. 
03. f = array([1,0,0,0,0,1,0,0,0])
04. print "f:",f
05. 
06. h = array([1,2,3])
07. print "h:",h
08. 
09. g1 = ia.iapconv2(f,h)
10. g2 = ia.iapconv2(h,f)
11. print "g1:",g1
12. print "g2:",g2
f: [1 0 0 0 0 1 0 0 0]
h: [1 2 3]
g1: [ 1.  2.  3.  0.  0.  1.  2.  3.  0.]
g2: [ 1.  2.  3.  0.  0.  1.  2.  3.  0.]

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.iapconv2(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):
[[ 1.  2.  3.  0.  0.  0.  0.  0.  0.]
 [ 4.  5.  6.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  1.  2.  3.  0.  0.  0.]
 [ 2.  3.  0.  4.  5.  6.  0.  0.  1.]
 [ 5.  6.  0.  0.  0.  0.  0.  0.  4.]]

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.iapconv2(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): 
[[[ 48.  42.  45.]
  [ 32.  26.  29.]
  [ 40.  34.  37.]]

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

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

Example with Image 2D

01. f = adread('cameraman.pgm')
02. adshow(f, title = 'a) - Original Image')
03. h = array([[-1,-1,-1],
04.            [ 0, 0, 0],
05.            [ 1, 1, 1]])
06. g = ia.iapconv2(f,h)
07. print "\nPrewitt´s Mask"
08. print h
09. 
10. gn = ia.ianormalize(g, [0,255])
11. adshow(gn, title = 'b) Prewitt´s Mask filtering')
12. 
13. 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.iapconv2(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

See also

  • iapconv - Periodic discrete convolution, kernel origin at center
  • iaconv - 2D or 3D linear discrete convolution.
  • iaptrans - Periodic translation.
  • iaconvteo - Illustrate the convolution theorem.

Contributions:

  • Francislei J. Silva (set 2013)
  • Roberto M. Souza (set 2013)