Function iainterpollin

Synopse

Perform linear interpolation

  • v = iainterpollin(f, pts)
    • f: 1D, 2D or 3D image to be interpolated
    • pts: array DxN with the points to interpolate (one in each column)
    • v: array of the interpolated values

For each cell the following equation needs to be solved:

1D: v = a + bx

2D: v = a + bx + cy + exy

3D: v = a + bx + cy + dz + exy + fxz + gyz + hxyz

where x, y and z are the fractional part of the coordinate.

001. def iainterpollin(f, pts):
002.     # f - one, two or three dimention array
003.     # pts - array of points to interpolate:
004.     # throws IndexError if there are indices out of range in pts
005.     from iainterpollin import iainterpollin1D, iainterpollin2D, iainterpollin3D
006.     import numpy as np
007.     f = f.astype(float)
008.     if f.ndim == 1:
009.         return iainterpollin1D(f, np.ravel(pts))
010.     if f.ndim == 2:
011.         return iainterpollin2D(f, pts)
012.     if f.ndim == 3:
013.         return iainterpollin3D(f, pts)
014. 
015. def iainterpollin1D(f, pts):
016.     # f - one dimention array
017.     # pts - array of points to interpolate:
018.     # throws IndexError if there are indices out of range in pts
019.     import numpy as np
020. 
021.     # integer indices
022.     ipts = np.floor(pts).astype(int)
023. 
024.     # fractional indices
025.     fpts = pts - ipts
026. 
027.     # workaround for the case some index equals the last valid index
028.     fpts[ipts>=f.shape[0]-1] += 1
029.     ipts[ipts>=f.shape[0]-1] -= 1
030. 
031.     I = ipts.copy()
032.     Ix = ipts.copy()+1
033. 
034.     a = f[I]
035.     b = f[Ix] - a
036. 
037.     return a + b*fpts
038. 
039. 
040. def iainterpollin2D(f, pts):
041.     # f - bi-dimentional image
042.     # pts - array (2xN) of points to interpolate:
043.     # throws IndexError if there are indices out of range in pts
044.     # let's consider x the first and y the second dimention
045.     import numpy as np
046. 
047.     # integer indices
048.     ipts = np.floor(pts).astype(int)
049. 
050.     # fractional indices
051.     fpts = pts - ipts
052. 
053.     # workaround for the case some index equals the last valid index in the respective dimension
054.     fpts[0][ipts[0]>=f.shape[0]-1] += 1
055.     ipts[0][ipts[0]>=f.shape[0]-1] -= 1
056.     fpts[1][ipts[1]>=f.shape[1]-1] += 1
057.     ipts[1][ipts[1]>=f.shape[1]-1] -= 1
058. 
059.     I = ipts.copy()
060.     Ix = ipts.copy()
061.     Ix[0] += 1
062.     Iy = ipts.copy()
063.     Iy[1] += 1
064.     Ixy = Ix.copy()
065.     Ixy[1] += 1
066. 
067.     a = f[I[0], I[1]]
068.     b = f[Ix[0], Ix[1]] - a
069.     c = f[Iy[0], Iy[1]] - a
070.     e = f[Ixy[0], Ixy[1]] - a - b - c
071. 
072.     return a + b*fpts[0] + c*fpts[1] + e*fpts[0]*fpts[1]
073. 
074. 
075. def iainterpollin3D(ff, pts):
076.     # ff - tri dimentional image
077.     # pts - array (3xN) of points to interpolate:
078.     # throws IndexError if there are indices out of range in pts
079.     # let's consider x the first dimension, y the second and z the third
080.     import numpy as np
081. 
082.     # integer indices
083.     ipts = np.floor(pts).astype(int)
084. 
085.     # fractional indices
086.     fpts = pts - ipts
087. 
088.     # workaround for the case some index equals the last valid index in the respective dimension
089.     fpts[0][ipts[0]>=ff.shape[0]-1] += 1
090.     ipts[0][ipts[0]>=ff.shape[0]-1] -= 1
091.     fpts[1][ipts[1]>=ff.shape[1]-1] += 1
092.     ipts[1][ipts[1]>=ff.shape[1]-1] -= 1
093.     fpts[2][ipts[2]>=ff.shape[2]-1] += 1
094.     ipts[2][ipts[2]>=ff.shape[2]-1] -= 1
095. 
096.     I = ipts.copy()
097.     Ix = ipts.copy()
098.     Ix[0] += 1
099.     Iy = ipts.copy()
100.     Iy[1] += 1
101.     Iz = ipts.copy()
102.     Iz[2] += 1
103.     Ixy = Ix.copy()
104.     Ixy[1] += 1
105.     Ixz = Ix.copy()
106.     Ixz[2] += 1
107.     Iyz = Iy.copy()
108.     Iyz[2] += 1
109.     Ixyz = Ixy.copy()
110.     Ixyz[2] += 1
111. 
112.     a = ff[I[0], I[1], I[2]]
113.     b = ff[Ix[0], Ix[1], Ix[2]] - a
114.     c = ff[Iy[0], Iy[1], Iy[2]] - a
115.     d = ff[Iz[0], Iz[1], Iz[2]] - a
116.     e = ff[Ixy[0], Ixy[1], Ixy[2]] - a - b - c
117.     f = ff[Ixz[0], Ixz[1], Ixz[2]] - a - b - d
118.     g = ff[Iyz[0], Iyz[1], Iyz[2]] - a - c - d
119.     h = ff[Ixyz[0], Ixyz[1], Ixyz[2]] - a - b - c - d - e - f - g
120. 
121.     return a + b*fpts[0] + c*fpts[1] + d*fpts[2] + e*fpts[0]*fpts[1] + f*fpts[0]*fpts[2] + g*fpts[1]*fpts[2] + h*fpts[0]*fpts[1]*fpts[2]

Example

01. import numpy as np
02. from ia636 import iainterpollin
03. 
04. print "1D interpolation"
05. Im = np.array([1, 7, 10, 5, 3, 0, -5, 6])
06. points = np.array([1.5, 2.1, 2.5, 2.9, 0.1])
07. print "Im= ", Im
08. print "points= ", points
09. print "Interpol vaules= ", iainterpollin(Im, points)
10. print "\n"
11. 
12. print "2D interpolation"
13. Im = np.arange(0, 9)
14. Im.shape = (3, 3)
15. print "Im=\n", Im
16. points = np.array([[0, 0.5, 0,   0.5, 1.5, 2],
17.                    [0, 0,   0.5, 0.5, 1.5, 2]])
18. print "points=\n", points
19. print "Interpol vaules= ", iainterpollin(Im, points)
20. print "\n"
21. 
22. print "3D interpolation"
23. Im = np.arange(0, 27)
24. Im.shape = (3, 3, 3)
25. print "Im=\n", Im
26. points = np.array([[0, 0.5, 0,   0, 0.5, 0.5, 2],
27.                    [0, 0,   0.5, 0, 0.5, 0.5, 2],
28.                    [0, 0,   0,   0.5, 0, 0.5, 2]])
29. print "points=\n", points
30. print "Interpol vaules= ", iainterpollin(Im, points)
1D interpolation
Im=  [ 1  7 10  5  3  0 -5  6]
points=  [ 1.5  2.1  2.5  2.9  0.1]
Interpol vaules=  [ 8.5  9.5  7.5  5.5  1.6]


2D interpolation
Im=
[[0 1 2]
 [3 4 5]
 [6 7 8]]
points=
[[ 0.   0.5  0.   0.5  1.5  2. ]
 [ 0.   0.   0.5  0.5  1.5  2. ]]
Interpol vaules=  [ 0.   1.5  0.5  2.   6.   8. ]


3D interpolation
Im=
[[[ 0  1  2]
  [ 3  4  5]
  [ 6  7  8]]

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

 [[18 19 20]
  [21 22 23]
  [24 25 26]]]
points=
[[ 0.   0.5  0.   0.   0.5  0.5  2. ]
 [ 0.   0.   0.5  0.   0.5  0.5  2. ]
 [ 0.   0.   0.   0.5  0.   0.5  2. ]]
Interpol vaules=  [  0.    4.5   1.5   0.5   6.    6.5  26. ]

See Also