Skip to content

Commit

Permalink
interpolation: changes in API (mostly backward compatible) and docs
Browse files Browse the repository at this point in the history 
in C++:
interpolate_value() renamed to trilinear_interpolation()
interpolate() renamed to interpolate_value()
default order=1 makes it backward compatible

the order argument: was 1/2/3, now it's 0/1/3
and corresponds to polynomial order.

in Python: the same, but interpolate() had no bindings
Also, added optional arg to_frac to grid.interpolate_position_array()

rewritten the Interpolation section in docs

it's a follow up to #324
functions that do Grid to Grid interpolation still need some work
  • Loading branch information
@wojdyr
wojdyr committed Sep 25, 2024
1 parent e98a6f4 commit 99cfabc
Show file tree
Hide file tree
Showing 4 changed files with 133 additions and 74 deletions.
130 changes: 88 additions & 42 deletions docs/grid.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,72 +342,100 @@ a bulk :ref:`solvent mask <solventmask>`.
Interpolation
-------------

To get a value corresponding to an arbitrary position,
you may use trilinear interpolation of the 8 nearest nodes,
or tricubic interpolation that uses 64 nodes.
Interpolation is used to obtain a value corresponding to an arbitrary position.
The most common interpolation methods are:

**C++**
* trilinear interpolation of the 8 nearest nodes,
* tricubic interpolation that uses 64 nodes.

::
Cubic interpolation is smoother than linear but may amplify noise.
This is illustrated in the plots below, which show density along two lines
in a grid filled with random numbers from [0, 1).
Trilinear interpolation is shown in blue, and tricubic -- in red.
The left plot shows density along a line in a random direction,
while the right plot shows density along a line parallel to one of the axes.

.. image:: img/trilinear-tricubic.png
:align: center
:scale: 100

T Grid<T>::interpolate_value(const Fractional& fctr) const
T Grid<T>::interpolate_value(const Position& ctr) const
In the functions below, the choice of interpolation is specified
using the `order` argument:

double Grid<T>::tricubic_interpolation(const Fractional& fctr) const
double Grid<T>::tricubic_interpolation(const Position& ctr) const
* order=0 -- the nearest grid point value,
* order=1 -- trilinear interpolation (default),
* order=3 -- tricubic interpolation.

// calculates also derivatives
**C++**

::

T Grid<T>::interpolate_value(const Fractional& fctr, int order=1) const
T Grid<T>::interpolate_value(const Position& ctr, int order=1) const
// You can also directly call the underlying functions
// trilinear_interpolation() and tricubic_interpolation().
// There is also a function that calculates derivatives:
std::array<double,4> Grid<T>::tricubic_interpolation_der(double x, double y, double z) const

**Python**

.. doctest::

>>> # trilinear interpolation
>>> grid.interpolate_value(gemmi.Fractional(1/24, 1/24, 1/24))
0.890625
>>> grid.interpolate_value(gemmi.Position(2, 3, 4))
2.0333263874053955
>>> grid.tricubic_interpolation(gemmi.Fractional(1/24, 1/24, 1/24))
>>> # nearest value
>>> grid.interpolate_value(gemmi.Fractional(1/24, 1/24, 1/24), order=0)
7.0
>>> grid.interpolate_value(gemmi.Position(2, 3, 4), order=0) # doctest: +ELLIPSIS
0.0
>>> # tricubic interpolation
>>> grid.interpolate_value(gemmi.Fractional(1/24, 1/24, 1/24), order=3)
1.283477783203125
>>> grid.tricubic_interpolation(gemmi.Position(2, 3, 4)) # doctest: +ELLIPSIS
2.6075661737715...
>>> # calculate also derivatives in directions of unit cell axes
>>> grid.interpolate_value(gemmi.Position(2, 3, 4), order=3) # doctest: +ELLIPSIS
2.607566...
>>> # calculate value and derivatives in the directions of unit cell axes
>>> grid.tricubic_interpolation_der(gemmi.Fractional(1/24, 1/24, 1/24))
[1.283477783203125, 35.523193359375, 36.343505859375, 35.523193359375]

The cubic interpolation is smoother than linear, but may amplify the noise.
This is illustrated on the plots below, which shows density along two lines
in a grid that was filled with random numbers from [0, 1).
Trilinear interpolation is blue, tricubic -- red.
The left plot shows density along a line in a random direction,
the right plot -- along a line parallel to one of the axes.

.. image:: img/trilinear-tricubic.png
:align: center
:scale: 100
*NumPy arrays*

*Implementation notes*
To interpolate the grid at positions listed as a NumPy array,
use `interpolate_position_array()`, which takes a NumPy array of positions
and returns a NumPy array of interpolated values:

Tricubic interpolation, as described
on `Wikipedia page <https://en.wikipedia.org/wiki/Tricubic_interpolation>`_ and in
`Appendix B of a PHENIX paper <https://journals.iucr.org/d/issues/2018/06/00/ic5103/#APPB>`_,
can be implemented either as 21 cubic interpolations, or using method
introduced by `Lekien & Marsen <https://doi.org/10.1002/nme.1296>`_ in 2005,
which involves 64x64 matrix of integral coefficients
(see also this `blog post <http://ianfaust.com/2016/03/20/Tricubic/>`_).
The latter method should be more efficient, but gemmi uses the former,
which takes ~100 ns. If you'd like to speed it up or to get derivatives,
contact developers.
.. doctest::
:skipif: numpy is None

>>> coords = numpy.array([[1, 2, 3], [2, 3, 4]], dtype=numpy.float32)
>>> grid.interpolate_position_array(coords)
array([0.4954875, 2.0333264], dtype=float32)

By default, it expects Cartesian coordinates and uses linear interpolation.
This can be changed by providing optional arguments `to_frac: Transform`
and `order: int`.

`to_frac` transforms input coordinates to fractional coordinates.
By default, it is a fractionalization matrix from the grid's unit cell.
If your coordinates are already fractional, pass the identity matrix:

.. doctest::
:skipif: numpy is None

*Optimization for Python*
>>> frac = numpy.array([[1/24, 1/24, 1/24]])
>>> grid.interpolate_position_array(frac, to_frac=gemmi.Transform())
array([0.890625], dtype=float32)

If you have a large number of points, making a Python function call
each time would be slow.
If these points are on a regular 3D grid (which may not be aligned
----

If the positions of interest are on a regular 3D grid (which may not be aligned
with our grid) call `interpolate_values()` (with s at the end)
with two arguments: a 3D NumPy array (for storing the results)
and a :ref:`Transform <transform>` that relates indices of the array
to positions in the grid:
and a :ref:`Transform <transform>` that relates the array's indices
to positions (in Angstroms) in the grid:

.. doctest::
:skipif: numpy is None
Expand All @@ -423,8 +451,26 @@ to positions in the grid:
>>> arr[10, 10, 10] # -> corresponds to Position(2, 3, 4)
2.0333264

(If your points are not on a regular grid -- get in touch -- there might be
another way.)
*One Grid to another*

For rescaling, rotating and translating maps, we have functions that use values
from one Grid to set values in another Grid.
Currently, the API is not ideal and should be revisited.
Therefore, for now we leave them undocumented:

gemmi::interpolate_grid()
gemmi::interpolate_grid_of_aligned_model2()

*Implementation note*

Tricubic interpolation, as described on the
`Wikipedia page <https://en.wikipedia.org/wiki/Tricubic_interpolation>`_ and in
`Appendix B of a PHENIX paper <https://journals.iucr.org/d/issues/2018/06/00/ic5103/#APPB>`_,
can be implemented either as 21 cubic interpolations or using a method
introduced by `Lekien & Marsden <https://doi.org/10.1002/nme.1296>`_ in 2005,
which involves 64x64 matrix of integral coefficients
(see also this `blog post <http://ianfaust.com/2016/03/20/Tricubic/>`_).
Gemmi uses the former method. It takes ~100 ns.

.. _masked_grid:

Expand Down
29 changes: 17 additions & 12 deletions include/gemmi/grid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -436,7 +436,7 @@ struct Grid : GridBase<T> {

/// https://en.wikipedia.org/wiki/Trilinear_interpolation
/// x,y,z are grid coordinates (x=1.5 is between 2nd and 3rd grid point).
T interpolate_value(double x, double y, double z) const {
T trilinear_interpolation(double x, double y, double z) const {
this->check_not_empty();
int u, v, w;
double xd = grid_modulo(x, nu, &u);
Expand All @@ -457,11 +457,11 @@ struct Grid : GridBase<T> {
}
return (T) lerp_(avg[0], avg[1], zd);
}
T interpolate_value(const Fractional& fctr) const {
return interpolate_value(fctr.x * nu, fctr.y * nv, fctr.z * nw);
T trilinear_interpolation(const Fractional& fctr) const {
return trilinear_interpolation(fctr.x * nu, fctr.y * nv, fctr.z * nw);
}
T interpolate_value(const Position& ctr) const {
return interpolate_value(unit_cell.fractionalize(ctr));
T trilinear_interpolation(const Position& ctr) const {
return trilinear_interpolation(unit_cell.fractionalize(ctr));
}

/// https://en.wikipedia.org/wiki/Tricubic_interpolation
Expand Down Expand Up @@ -541,14 +541,17 @@ struct Grid : GridBase<T> {
copy[i][j][k] = this->get_value_q(u_indices[i], v_indices[j], w_indices[k]);
}

/// @param order 1=nearest, 2=linear, 3=cubic interpolation
T interpolate(const Fractional& f, int order) const {
/// @param order 0=nearest, 1=linear, 3=cubic interpolation
T interpolate_value(const Fractional& f, int order=1) const {
switch (order) {
case 1: return *const_cast<Grid<T>*>(this)->get_nearest_point(f).value;
case 2: return interpolate_value(f);
case 0: return *const_cast<Grid<T>*>(this)->get_nearest_point(f).value;
case 1: return trilinear_interpolation(f);
case 3: return (T) tricubic_interpolation(f);
}
throw std::invalid_argument("interpolation \"order\" must 1, 2 or 3");
throw std::invalid_argument("interpolation \"order\" must 0, 1 or 3");
}
T interpolate_value(const Position& ctr, int order=1) const {
return interpolate_value(unit_cell.fractionalize(ctr), order);
}

void get_subarray(T* dest, std::array<int,3> start, std::array<int,3> shape) const {
Expand Down Expand Up @@ -786,16 +789,18 @@ struct Grid : GridBase<T> {
};

// TODO: add argument Box<Fractional> src_extent
// cf. interpolate_grid_of_aligned_model2() in solmask.hpp
// cf interpolate_values in python/grid.cpp
template<typename T>
void interpolate_grid(Grid<T>& dest, const Grid<T>& src, const Transform& tr, int order=2) {
void interpolate_grid(Grid<T>& dest, const Grid<T>& src, const Transform& tr, int order=1) {
FTransform frac_tr = src.unit_cell.frac.combine(tr).combine(dest.unit_cell.orth);
size_t idx = 0;
for (int w = 0; w != dest.nw; ++w)
for (int v = 0; v != dest.nv; ++v)
for (int u = 0; u != dest.nu; ++u, ++idx) {
Fractional dest_fr = dest.get_fractional(u, v, w);
Fractional src_fr = frac_tr.apply(dest_fr);
dest.data[idx] = src.interpolate(src_fr, order);
dest.data[idx] = src.interpolate_value(src_fr, order);
}
}

Expand Down
4 changes: 2 additions & 2 deletions include/gemmi/solmask.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -396,7 +396,7 @@ template<typename T>
void interpolate_grid_of_aligned_model2(Grid<T>& dest, const Grid<T>& src,
const Transform& tr,
const Model& dest_model, double radius,
int order=2) {
int order=1) {
Grid<NodeInfo> mask;
mask.copy_metadata_from(dest);
mask_with_node_info(mask, dest_model, radius);
Expand All @@ -408,7 +408,7 @@ void interpolate_grid_of_aligned_model2(Grid<T>& dest, const Grid<T>& src,
if (ni.found) {
Fractional dest_fr = dest.get_fractional(ni.u, ni.v, ni.w);
Fractional src_fr = frac_tr.apply(dest_fr);
dest.data[idx] = src.interpolate(src_fr, order);
dest.data[idx] = src.interpolate_value(src_fr, order);
}
}
}
Expand Down
44 changes: 26 additions & 18 deletions python/grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,23 +162,38 @@ void add_grid_interpolation(nb::class_<Grid<T>, GridBase<T>>& grid) {
using Gr = Grid<T>;
grid
.def("interpolate_value",
(T (Gr::*)(const Fractional&) const) &Gr::interpolate_value)
(T (Gr::*)(const Fractional&, int) const) &Gr::interpolate_value,
nb::arg(), nb::arg("order")=1)
.def("interpolate_value",
(T (Gr::*)(const Position&) const) &Gr::interpolate_value)
(T (Gr::*)(const Position&, int) const) &Gr::interpolate_value,
nb::arg(), nb::arg("order")=1)
// deprecated, use interpolate_value(..., order=3)
.def("tricubic_interpolation",
(double (Gr::*)(const Fractional&) const) &Gr::tricubic_interpolation)
// deprecated, use interpolate_value(..., order=3)
.def("tricubic_interpolation",
(double (Gr::*)(const Position&) const) &Gr::tricubic_interpolation)
.def("tricubic_interpolation_der",
(std::array<double,4> (Gr::*)(const Fractional&) const)
&Gr::tricubic_interpolation_der)
.def("interpolate_position_array",
[](const Gr& self, nb::ndarray<T, nb::shape<-1,3>, nb::device::cpu> xyz, int order) {
[](const Gr& self, nb::ndarray<double, nb::shape<-1,3>, nb::device::cpu> xyz,
int order, const Transform* to_frac) {
auto xyz_view = xyz.view();
size_t len = xyz_view.shape(0);
auto values = make_numpy_array<T>({len});
T* data = values.data();
const Transform& frac = to_frac ? *to_frac : self.unit_cell.frac;
for (size_t i = 0; i < len; ++i) {
Position pos(xyz_view(i, 0), xyz_view(i, 1), xyz_view(i, 2));
Fractional fpos = self.unit_cell.fractionalize(pos);
data[i] = self.interpolate(fpos, order);
Fractional fpos = Fractional(frac.apply(pos));
data[i] = self.interpolate_value(fpos, order);
}
return values;
}, nb::arg(), nb::arg("order")=2)
// TODO: find a better name for this func, perhaps interpolate_array?
}, nb::arg("xyz"), nb::arg("order")=1, nb::arg("to_frac")=nb::none())
// The name of this function is not very descriptive, but since it's used
// in a few external projects, renaming it isn't worth the hassle.
// cf. interpolate_grid
.def("interpolate_values",
[](const Gr& self, nb::ndarray<nb::numpy, T, nb::ndim<3>> arr,
const Transform& tr, int order) {
Expand All @@ -188,16 +203,9 @@ void add_grid_interpolation(nb::class_<Grid<T>, GridBase<T>>& grid) {
for (size_t k = 0; k < r.shape(2); ++k) {
Position pos(tr.apply(Vec3(i, j, k)));
Fractional fpos = self.unit_cell.fractionalize(pos);
r(i, j, k) = self.interpolate(fpos, order);
r(i, j, k) = self.interpolate_value(fpos, order);
}
}, nb::arg().noconvert(), nb::arg(), nb::arg("order")=2)
.def("tricubic_interpolation",
(double (Gr::*)(const Fractional&) const) &Gr::tricubic_interpolation)
.def("tricubic_interpolation",
(double (Gr::*)(const Position&) const) &Gr::tricubic_interpolation)
.def("tricubic_interpolation_der",
(std::array<double,4> (Gr::*)(const Fractional&) const)
&Gr::tricubic_interpolation_der)
}, nb::arg().noconvert(), nb::arg(), nb::arg("order")=1)
;
}

Expand Down Expand Up @@ -265,10 +273,10 @@ void add_grid(nb::module_& m) {
.def("set_to_zero", &SolventMasker::set_to_zero)
;
m.def("interpolate_grid", &interpolate_grid<float>,
nb::arg("dest"), nb::arg("src"), nb::arg("tr"), nb::arg("order")=2);
nb::arg("dest"), nb::arg("src"), nb::arg("tr"), nb::arg("order")=1);
m.def("interpolate_grid_of_aligned_model2", &interpolate_grid_of_aligned_model2<float>,
nb::arg("dest"), nb::arg("src"), nb::arg("tr"),
nb::arg("dest_model"), nb::arg("radius"), nb::arg("order")=2);
nb::arg("dest_model"), nb::arg("radius"), nb::arg("order")=1);


// from blob.hpp
Expand Down

0 comments on commit 99cfabc

Please sign in to comment.