[theano-users] Interpolating a 2D array using theano variables.

Discussion:

Vital Fernández

2018-05-23 16:20:16 UTC

I recently started using theano within the pymc3 library.

I wonder if anyone could point me towards the right direction in this issue:

I have a 2D array (lets say a map with X,Y coordinantes and the matrix
values represents the altitude).

For a given set of coordinates (theano variables) I want to get the scalar
interpolated from the 4 nearest points in the grid. (in which X and Y are
ordered)

In numpy there are many ways to do this but I have not found and
equivalent function in theano. Should I get the indeces manually (for
example using numpy functions: idx_X = np.argmin(np.abs(X_array - x_cord)))

Thank you for any advice

--
---
You received this message because you are subscribed to the Google Groups "theano-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to theano-users+***@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Robb Brown

2018-05-23 19:48:08 UTC

Permalink

It's fairly easy to write an interpolation routine in Theano.

Here are nearest neighbour (3d) and trilinear for example. For your 2D
case you can just knock off a variable. "image" is a tensor containing the
image data, points is a list of coordinates giving the points you'd like
the interpolated value of. I expect these will work, although be warned, I
haven't tested them since Theano 0.7 or 0.8.

Nearest neighbour:

x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)

nearestNeighbour =
image[T.cast(T.floor(z),'uint16'),T.cast(T.floor(y),'uint16'),T.cast(T.floor(x),'uint16')]
nearestNeighbour.name = self.name + '-nearestNeighbourSampler'

Trilinear:

x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)

x0 = T.cast(T.floor(x),'uint16')
y0 = T.cast(T.floor(y),'uint16')
z0 = T.cast(T.floor(z),'uint16')

x1 = x0+1
y1 = y0+1
z1 = z0+1

xd = (x-x0)
yd = (y-y0)
zd = (z-z0)

c00 = image[z0,y0,x0]*(1-xd) + image[z0,y0,x1]*xd
c10 = image[z0,y1,x0]*(1-xd) + image[z0,y1,x1]*xd
c01 = image[z1,y0,x0]*(1-xd) + image[z1,y0,x1]*xd
c11 = image[z1,y1,x0]*(1-xd) + image[z1,y1,x1]*xd

c0 = c00*(1-yd) + c10*yd
c1 = c01*(1-yd) + c11*yd

trilinear = c0*(1-zd) + c1*zd

Vital Fernández

2018-05-23 20:00:12 UTC

Permalink

Wow Robb. Thank you very much. I am going to play with this to make a
complete script then I will post it here :).

Thanks

Post by Robb Brown
It's fairly easy to write an interpolation routine in Theano.
Here are nearest neighbour (3d) and trilinear for example. For your 2D
case you can just knock off a variable. "image" is a tensor containing the
image data, points is a list of coordinates giving the points you'd like
the interpolated value of. I expect these will work, although be warned, I
haven't tested them since Theano 0.7 or 0.8.
x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)
nearestNeighbour =
image[T.cast(T.floor(z),'uint16'),T.cast(T.floor(y),'uint16'),T.cast(T.floor(x),'uint16')]
nearestNeighbour.name = self.name + '-nearestNeighbourSampler'
x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)
x0 = T.cast(T.floor(x),'uint16')
y0 = T.cast(T.floor(y),'uint16')
z0 = T.cast(T.floor(z),'uint16')
x1 = x0+1
y1 = y0+1
z1 = z0+1
xd = (x-x0)
yd = (y-y0)
zd = (z-z0)
c00 = image[z0,y0,x0]*(1-xd) + image[z0,y0,x1]*xd
c10 = image[z0,y1,x0]*(1-xd) + image[z0,y1,x1]*xd
c01 = image[z1,y0,x0]*(1-xd) + image[z1,y0,x1]*xd
c11 = image[z1,y1,x0]*(1-xd) + image[z1,y1,x1]*xd
c0 = c00*(1-yd) + c10*yd
c1 = c01*(1-yd) + c11*yd
trilinear = c0*(1-zd) + c1*zd

Robb Brown

2018-05-23 20:26:58 UTC

Permalink

Youâre very welcome.

Robert Brown, PhD
McConnell Brain Imaging Centre
Montreal Neurological Institute
Montreal, Canada
www.robbtech.com <http://www.robbtech.com/>
github.com/robb-brown <http://github.com/robb-brown>

Wow Robb. Thank you very much. I am going to play with this to make a complete script then I will post it here :).
Thanks
It's fairly easy to write an interpolation routine in Theano.
Here are nearest neighbour (3d) and trilinear for example. For your 2D case you can just knock off a variable. "image" is a tensor containing the image data, points is a list of coordinates giving the points you'd like the interpolated value of. I expect these will work, although be warned, I haven't tested them since Theano 0.7 or 0.8.
x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)
nearestNeighbour = image[T.cast(T.floor(z),'uint16'),T.cast(T.floor(y),'uint16'),T.cast(T.floor(x),'uint16')]
nearestNeighbour.name = self.name <http://self.name/> + '-nearestNeighbourSampler'
x = T.clip(points[2],0,T.shape(image)[2]-2)
y = T.clip(points[1],0,T.shape(image)[1]-2)
z = T.clip(points[0],0,T.shape(image)[0]-2)
x0 = T.cast(T.floor(x),'uint16')
y0 = T.cast(T.floor(y),'uint16')
z0 = T.cast(T.floor(z),'uint16')
x1 = x0+1
y1 = y0+1
z1 = z0+1
xd = (x-x0)
yd = (y-y0)
zd = (z-z0)
c00 = image[z0,y0,x0]*(1-xd) + image[z0,y0,x1]*xd
c10 = image[z0,y1,x0]*(1-xd) + image[z0,y1,x1]*xd
c01 = image[z1,y0,x0]*(1-xd) + image[z1,y0,x1]*xd
c11 = image[z1,y1,x0]*(1-xd) + image[z1,y1,x1]*xd
c0 = c00*(1-yd) + c10*yd
c1 = c01*(1-yd) + c11*yd
trilinear = c0*(1-zd) + c1*zd
--
---
You received this message because you are subscribed to a topic in the Google Groups "theano-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/theano-users/lPN_vWkOIK0/unsubscribe <https://groups.google.com/d/topic/theano-users/lPN_vWkOIK0/unsubscribe>.
For more options, visit https://groups.google.com/d/optout <https://groups.google.com/d/optout>.

Nicolas Budyn

2018-07-12 11:20:35 UTC

Permalink

I found this thread very useful, thank you both. It took me time to get the
code working with theano 1.0.1 and numpy 1.14.3 because I was obtaining the
error "ValueError: expected an ndarray" when evaluating the trilinear
interpolator. The issue was coming from the "uint16" type used for
indexing: using "int32" instead, it works flawlessly. Also for the nearest
interpolation, "uint16" produces an AdvancedSubtensor node which in this
case around 100 times slower than the Subtensor node produced when using
"int32". In conclusion, I used "int32" everywhere. The complete code is
below.

Regards

Nicolas

#########
import theano
import theano.tensor as T

int_ = "int32"

# Init
image = T.dtensor3("image")
points = T.dvector("points")

x = T.clip(points[2], 0, T.shape(image)[2] - 2)
y = T.clip(points[1], 0, T.shape(image)[1] - 2)
z = T.clip(points[0], 0, T.shape(image)[0] - 2)

# Nearest neighbour:
nearest = image[
T.cast(T.floor(z), int_), T.cast(T.floor(y), int_), T.cast(T.floor(x),
int_)
]
nearest_f = theano.function([image, points], nearest)

# Trilinear:
x0 = T.cast(T.floor(x), int_)
y0 = T.cast(T.floor(y), int_)
z0 = T.cast(T.floor(z), int_)

x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1

xd = x - x0
yd = y - y0
zd = z - z0

c00 = image[z0, y0, x0] * (1 - xd) + image[z0, y0, x1] * xd
c10 = image[z0, y1, x0] * (1 - xd) + image[z0, y1, x1] * xd
c01 = image[z1, y0, x0] * (1 - xd) + image[z1, y0, x1] * xd
c11 = image[z1, y1, x0] * (1 - xd) + image[z1, y1, x1] * xd

c0 = c00 * (1 - yd) + c10 * yd
c1 = c01 * (1 - yd) + c11 * yd

trilinear = c0 * (1 - zd) + c1 * zd
trilinear_f = theano.function([image, points], trilinear)

# Test
import numpy as np

zz, yy, xx = np.mgrid[0:10, 0:10, 0:10]
zz = zz.astype(np.float_)
yy = yy.astype(np.float_)
xx = xx.astype(np.float_)
image_val = xx ** 2 + yy
image_val = image_val.astype(np.float_)
points_val = np.array([0., 6., 9.])

print(nearest_f(image_val, points_val))
print(trilinear_f(image_val, points_val))
#########

Post by Robb Brown
Youâre very welcome.
Robert Brown, PhD
McConnell Brain Imaging Centre
Montreal Neurological Institute
Montreal, Canada
www.robbtech.com
github.com/robb-brown
Wow Robb. Thank you very much. I am going to play with this to make a
complete script then I will post it here :).
Thanks

--
---
You received this message because you are subscribed to a topic in the
Google Groups "theano-users" group.
To unsubscribe from this topic, visit
https://groups.google.com/d/topic/theano-users/lPN_vWkOIK0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to
For more options, visit https://groups.google.com/d/optout.

Vital Fernández

2018-07-12 14:50:02 UTC

Permalink

Thank you Nicolas very much for sharing it is very neat.

Actually you can actually see on the forums and on github posts about
theano having an ops for bilinear interpolation but I have not seen any
working example... I think I should post a new topic to check

Post by Nicolas Budyn
I found this thread very useful, thank you both. It took me time to get
the code working with theano 1.0.1 and numpy 1.14.3 because I was obtaining
the error "ValueError: expected an ndarray" when evaluating the trilinear
interpolator. The issue was coming from the "uint16" type used for
indexing: using "int32" instead, it works flawlessly. Also for the nearest
interpolation, "uint16" produces an AdvancedSubtensor node which in this
case around 100 times slower than the Subtensor node produced when using
"int32". In conclusion, I used "int32" everywhere. The complete code is
below.
Regards
Nicolas
#########
import theano
import theano.tensor as T
int_ = "int32"
# Init
image = T.dtensor3("image")
points = T.dvector("points")
x = T.clip(points[2], 0, T.shape(image)[2] - 2)
y = T.clip(points[1], 0, T.shape(image)[1] - 2)
z = T.clip(points[0], 0, T.shape(image)[0] - 2)
nearest = image[
T.cast(T.floor(z), int_), T.cast(T.floor(y), int_), T.cast(T.floor(x),
int_)
]
nearest_f = theano.function([image, points], nearest)
x0 = T.cast(T.floor(x), int_)
y0 = T.cast(T.floor(y), int_)
z0 = T.cast(T.floor(z), int_)
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
xd = x - x0
yd = y - y0
zd = z - z0
c00 = image[z0, y0, x0] * (1 - xd) + image[z0, y0, x1] * xd
c10 = image[z0, y1, x0] * (1 - xd) + image[z0, y1, x1] * xd
c01 = image[z1, y0, x0] * (1 - xd) + image[z1, y0, x1] * xd
c11 = image[z1, y1, x0] * (1 - xd) + image[z1, y1, x1] * xd
c0 = c00 * (1 - yd) + c10 * yd
c1 = c01 * (1 - yd) + c11 * yd
trilinear = c0 * (1 - zd) + c1 * zd
trilinear_f = theano.function([image, points], trilinear)
# Test
import numpy as np
zz, yy, xx = np.mgrid[0:10, 0:10, 0:10]
zz = zz.astype(np.float_)
yy = yy.astype(np.float_)
xx = xx.astype(np.float_)
image_val = xx ** 2 + yy
image_val = image_val.astype(np.float_)
points_val = np.array([0., 6., 9.])
print(nearest_f(image_val, points_val))
print(trilinear_f(image_val, points_val))
#########

Robb Brown

2018-07-13 22:29:59 UTC

Permalink

Glad you found it useful. Thanks for the tip about the data type.

Robert Brown, PhD
McConnell Brain Imaging Centre
Montreal Neurological Institute
Montreal, Canada
www.robbtech.com <http://www.robbtech.com/>
github.com/robb-brown <http://github.com/robb-brown>

I found this thread very useful, thank you both. It took me time to get the code working with theano 1.0.1 and numpy 1.14.3 because I was obtaining the error "ValueError: expected an ndarray" when evaluating the trilinear interpolator. The issue was coming from the "uint16" type used for indexing: using "int32" instead, it works flawlessly. Also for the nearest interpolation, "uint16" produces an AdvancedSubtensor node which in this case around 100 times slower than the Subtensor node produced when using "int32". In conclusion, I used "int32" everywhere. The complete code is below.
Regards
Nicolas
#########
import theano
import theano.tensor as T
int_ = "int32"
# Init
image = T.dtensor3("image")
points = T.dvector("points")
x = T.clip(points[2], 0, T.shape(image)[2] - 2)
y = T.clip(points[1], 0, T.shape(image)[1] - 2)
z = T.clip(points[0], 0, T.shape(image)[0] - 2)
nearest = image[
T.cast(T.floor(z), int_), T.cast(T.floor(y), int_), T.cast(T.floor(x), int_)
]
nearest_f = theano.function([image, points], nearest)
x0 = T.cast(T.floor(x), int_)
y0 = T.cast(T.floor(y), int_)
z0 = T.cast(T.floor(z), int_)
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
xd = x - x0
yd = y - y0
zd = z - z0
c00 = image[z0, y0, x0] * (1 - xd) + image[z0, y0, x1] * xd
c10 = image[z0, y1, x0] * (1 - xd) + image[z0, y1, x1] * xd
c01 = image[z1, y0, x0] * (1 - xd) + image[z1, y0, x1] * xd
c11 = image[z1, y1, x0] * (1 - xd) + image[z1, y1, x1] * xd
c0 = c00 * (1 - yd) + c10 * yd
c1 = c01 * (1 - yd) + c11 * yd
trilinear = c0 * (1 - zd) + c1 * zd
trilinear_f = theano.function([image, points], trilinear)
# Test
import numpy as np
zz, yy, xx = np.mgrid[0:10, 0:10, 0:10]
zz = zz.astype(np.float_)
yy = yy.astype(np.float_)
xx = xx.astype(np.float_)
image_val = xx ** 2 + yy
image_val = image_val.astype(np.float_)
points_val = np.array([0., 6., 9.])
print(nearest_f(image_val, points_val))
print(trilinear_f(image_val, points_val))
#########
Youâre very welcome.
Robert Brown, PhD
McConnell Brain Imaging Centre
Montreal Neurological Institute
Montreal, Canada
www.robbtech.com <http://www.robbtech.com/>
github.com/robb-brown <http://github.com/robb-brown>

--
---
You received this message because you are subscribed to a topic in the Google Groups "theano-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/theano-users/lPN_vWkOIK0/unsubscribe <https://groups.google.com/d/topic/theano-users/lPN_vWkOIK0/unsubscribe>.
For more options, visit https://groups.google.com/d/optout <https://groups.google.com/d/optout>.