From 9c9a8db46ff6659578fe024d9e252efdcf5f1092 Mon Sep 17 00:00:00 2001 From: Nicola Demo Date: Mon, 15 May 2017 15:35:01 +0200 Subject: [PATCH] Add IDW deformation technique - add the IDW class to handle the deformation - add the IDWParameters class - add test about parameters and deformation - some minor changes in rbf deformation --- docs/source/code.rst | 3 +- docs/source/idw.rst | 20 +++ docs/source/params_idw.rst | 21 +++ pygem/__init__.py | 20 +-- pygem/idw.py | 121 ++++++++++++++++ pygem/params_idw.py | 131 ++++++++++++++++++ pygem/radial.py | 17 ++- .../test_datasets/parameters_idw_default.prm | 34 +++++ tests/test_datasets/parameters_idw_deform.prm | 34 +++++ tests/test_idw.py | 38 +++++ tests/test_idwparams.py | 95 +++++++++++++ 11 files changed, 516 insertions(+), 18 deletions(-) create mode 100644 docs/source/idw.rst create mode 100644 docs/source/params_idw.rst create mode 100644 pygem/idw.py create mode 100644 pygem/params_idw.py create mode 100644 tests/test_datasets/parameters_idw_default.prm create mode 100644 tests/test_datasets/parameters_idw_deform.prm create mode 100644 tests/test_idw.py create mode 100644 tests/test_idwparams.py diff --git a/docs/source/code.rst b/docs/source/code.rst index a09ad4c..8b0310c 100644 --- a/docs/source/code.rst +++ b/docs/source/code.rst @@ -20,4 +20,5 @@ Code Documentation stephandler utils gui - + idw + params_idw diff --git a/docs/source/idw.rst b/docs/source/idw.rst new file mode 100644 index 0000000..08b89c6 --- /dev/null +++ b/docs/source/idw.rst @@ -0,0 +1,20 @@ +Inverse Distance Weighting +========================== + +.. currentmodule:: pygem.idw + +.. automodule:: pygem.idw + +.. autosummary:: + :toctree: _summaries + :nosignatures: + + IDW.perform + +.. autoclass:: IDW + :members: + :private-members: + :undoc-members: + :show-inheritance: + :noindex: + diff --git a/docs/source/params_idw.rst b/docs/source/params_idw.rst new file mode 100644 index 0000000..489b34e --- /dev/null +++ b/docs/source/params_idw.rst @@ -0,0 +1,21 @@ +IDWParameters +================= + +.. currentmodule:: pygem.params_idw + +.. automodule:: pygem.params_idw + +.. autosummary:: + :toctree: _summaries + :nosignatures: + + IDWParameters.__str__ + IDWParameters.read_parameters + IDWParameters.write_parameters + +.. autoclass:: IDWParameters + :members: + :private-members: + :undoc-members: + :show-inheritance: + :noindex: diff --git a/pygem/__init__.py b/pygem/__init__.py index 5ec0b1a..43ae037 100644 --- a/pygem/__init__.py +++ b/pygem/__init__.py @@ -1,21 +1,25 @@ +""" +""" __all__ = [ 'affine', 'filehandler', 'freeform', 'radial', 'openfhandler', 'params', 'stlhandler', 'unvhandler', 'vtkhandler', 'nurbshandler', 'stephandler', - 'igeshandler', 'utils', 'gui', 'khandler' + 'igeshandler', 'utils', 'gui', 'khandler', 'idw', 'params_idw' ] from . import affine from . import filehandler from . import freeform -from . import gui -from . import igeshandler -from . import khandler -from . import nurbshandler +from . import radial from . import openfhandler from . import params -from . import radial -from . import stephandler from . import stlhandler from . import unvhandler -from . import utils from . import vtkhandler +from . import nurbshandler +from . import stephandler +from . import igeshandler +from . import utils +from . import gui +from . import khandler +from . import idw +from . import params_idw diff --git a/pygem/idw.py b/pygem/idw.py new file mode 100644 index 0000000..7293cb8 --- /dev/null +++ b/pygem/idw.py @@ -0,0 +1,121 @@ +""" +Module focused on the Inverse Distance Weighting interpolation technique. +The IDW algorithm is an average moving interpolation that is usually applied to +highly variable data. The main idea of this interpolation strategy lies in +fact that it is not desirable to honour local high/low values but rather to look +at a moving average of nearby data points and estimate the local trends. +The node value is calculated by averaging the weighted sum of all the points. +Data points that lie progressively farther from the node inuence much less the +computed value than those lying closer to the node. + +:Theoretical Insight: + + This implementation is based on the simplest form of inverse distance + weighting interpolation, proposed by D. Shepard, A two-dimensional + interpolation function for irregularly-spaced data, Proceedings of the 23 rd + ACM National Conference. + + The interpolation value :math:`u` of a given point :math:`\\mathrm{x}` + from a set of samples :math:`u_k = u(\\mathrm{x}_k)`, with + :math:`k = 1,2,\dotsc,\\mathcal{N}`, is given by: + + .. math:: + u(\\mathrm{x}) = \\displaystyle\\sum_{k=1}^\\mathcal{N} + \\frac{w(\\mathrm{x},\\mathrm{x}_k)} + {\\displaystyle\\sum_{j=1}^\\mathcal{N} w(\\mathrm{x},\\mathrm{x}_j)} + u_k + + + where, in general, :math:`w(\\mathrm{x}, \\mathrm{x}_i)` represents the + weighting function: + + .. math:: + w(\\mathrm{x}, \\mathrm{x}_i) = \\| \\mathrm{x} - \\mathrm{x}_i \\|^{-p} + + being :math:`\\| \\mathrm{x} - \\mathrm{x}_i \\|^{-p} \\ge 0` is the + Euclidean distance between :math:`\\mathrm{x}` and data point + :math:`\\mathrm{x}_i` and :math:`p` is a power parameter, typically equal to + 2. + +""" +import numpy as np + +from scipy.spatial.distance import cdist + + +class IDW(object): + """ + Class that handles the IDW technique. + + :param idw_parameters: the parameters of the IDW + :type idw_parameters: :class:`IDWParameters` + :param numpy.ndarray original_mesh_points: coordinates of the original + points of the mesh. + + :cvar parameters: the parameters of the IDW. + :vartype parameters: :class:`~pygem.params_idw.IDWParameters` + :cvar numpy.ndarray original_mesh_points: coordinates of the original + points of the mesh. + :cvar numpy.ndarray modified_mesh_points: coordinates of the deformed + points of the mesh. + + :Example: + + >>> from pygem.idw import IDW + >>> from pygem.params_idw import IDWParameters + >>> import numpy as np + >>> params = IDWParameters() + >>> params.read_parameters('tests/test_datasets/parameters_idw_cube.prm') + >>> nx, ny, nz = (20, 20, 20) + >>> mesh = np.zeros((nx * ny * nz, 3)) + >>> xv = np.linspace(0, 1, nx) + >>> yv = np.linspace(0, 1, ny) + >>> zv = np.linspace(0, 1, nz) + >>> z, y, x = np.meshgrid(zv, yv, xv) + >>> mesh = np.array([x.ravel(), y.ravel(), z.ravel()]) + >>> original_mesh_points = mesh.T + >>> idw = IDW(rbf_parameters, original_mesh_points) + >>> idw.perform() + >>> new_mesh_points = idw.modified_mesh_points + """ + + def __init__(self, idw_parameters, original_mesh_points): + self.parameters = idw_parameters + self.original_mesh_points = original_mesh_points + self.modified_mesh_points = None + + def perform(self): + """ + This method performs the deformation of the mesh points. After the + execution it sets `self.modified_mesh_points`. + """ + + def distance(u, v): + return np.linalg.norm(u - v, self.parameters.power) + + # Compute displacement of the control points + displ = ( + self.parameters.deformed_control_points - + self.parameters.original_control_points + ) + + # Compute the distance between the mesh points and the control points + dist = cdist( + self.original_mesh_points, + self.parameters.original_control_points, + distance + ) + + # Weights are set as the reciprocal of the distance if the distance is + # not zero, otherwise 1.0 where distance is zero. + weights = np.zeros(dist.shape) + for i, d in enumerate(dist): + weights[i] = 1. / d if d.all() else np.where(d == 0.0, 1.0, 0.0) + + + offset = np.array([ + np.sum(displ * wi[:, np.newaxis] / np.sum(wi), axis=0) + for wi in weights + ]) + + self.modified_mesh_points = self.original_mesh_points + offset diff --git a/pygem/params_idw.py b/pygem/params_idw.py new file mode 100644 index 0000000..8dd0fe6 --- /dev/null +++ b/pygem/params_idw.py @@ -0,0 +1,131 @@ +import numpy as np +import os +try: + import configparser as configparser +except ImportError: + import ConfigParser as configparser + +class IDWParameters(object): + """ + Class that handles the Inverse Distance Weighting parameters in terms of + control points. + + :cvar int p: the power parameter. The default value is 2. + :cvar numpy.ndarray original_control_points: it is an + `n_control_points`-by-3 array with the coordinates of the original + interpolation control points before the deformation. The default is the + vertices of the unit cube. + :cvar numpy.ndarray deformed_control_points: it is an + `n_control_points`-by-3 array with the coordinates of the interpolation + control points after the deformation. The default is the vertices of + the unit cube. + """ + def __init__(self): + self.power = 2 + self.original_control_points = np.array([ + [0., 0., 0.], + [0., 0., 1.], + [0., 1., 0.], + [1., 0., 0.], + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], + [1., 1., 1.] + ]) + self.deformed_control_points = np.array([ + [0., 0., 0.], + [0., 0., 1.], + [0., 1., 0.], + [1., 0., 0.], + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], + [1., 1., 1.] + ]) + + + def read_parameters(self, filename): + """ + Reads in the parameters file and fill the self structure. + + :param string filename: parameters file to be read in. + """ + if not isinstance(filename, str): + raise TypeError('filename must be a string') + + if not os.path.isfile(filename): + raise IOError('filename does not exist') + + config = configparser.RawConfigParser() + config.read(filename) + + self.power = config.getint('Inverse Distance Weighting', 'power') + + ctrl_points = config.get('Control points', 'original control points') + self.original_control_points = np.array( + [line.split() for line in ctrl_points.split('\n')], + dtype=float + ) + + defo_points = config.get('Control points', 'deformed control points') + self.deformed_control_points = np.array( + [line.split() for line in defo_points.split('\n')], + dtype=float + ) + + + def write_parameters(self, filename): + """ + This method writes a parameters file (.prm) called `filename` and fills + it with all the parameters class members. + + :param string filename: parameters file to be written out. + """ + if not isinstance(filename, str): + raise TypeError("filename must be a string") + + output_string = "" + output_string += "\n[Inverse Distance Weighting]\n" + output_string += "# This section describes the settings of idw.\n\n" + output_string += "# the power parameter\n" + output_string += "power = {}\n".format(self.power) + + output_string += "\n\n[Control points]\n" + output_string += "# This section describes the IDW control points.\n\n" + output_string += "# original control points collects the coordinates\n" + output_string += "# of the interpolation control points before the\n" + output_string += "# deformation.\n" + + output_string += "original control points: " + output_string += ( + ' '.join(map(str, self.original_control_points[0])) + "\n" + ) + for points in self.original_control_points[1:]: + output_string += 25 * ' ' + ' '.join(map(str, points)) + "\n" + output_string += "\n" + output_string += "# deformed control points collects the coordinates\n" + output_string += "# of the interpolation control points after the\n" + output_string += "# deformation.\n" + output_string += "deformed control points: " + output_string += ( + ' '.join(map(str, self.original_control_points[0])) + "\n" + ) + for points in self.deformed_control_points[1:]: + output_string += 25 * ' ' + ' '.join(map(str, points)) + "\n" + + with open(filename, 'w') as f: + f.write(output_string) + + + def __str__(self): + """ + This method prints all the IDW parameters on the screen. Its purpose is + for debugging. + """ + string = '' + string += 'p = {}\n'.format(self.power) + string += '\noriginal_control_points =\n' + string += '{}\n'.format(self.original_control_points) + string += '\ndeformed_control_points =\n' + string += '{}\n'.format(self.deformed_control_points) + return string diff --git a/pygem/radial.py b/pygem/radial.py index b389565..ab48a84 100644 --- a/pygem/radial.py +++ b/pygem/radial.py @@ -47,6 +47,8 @@ """ import numpy as np +from scipy.spatial.distance import cdist + class RBF(object): """ @@ -223,11 +225,9 @@ def _distance_matrix(self, X1, X2): :return: matrix: the matrix D. :rtype: numpy.ndarray """ - m, n = X1.shape[0], X2.shape[0] - matrix = np.zeros(shape=(m, n)) - for i in range(0, m): - for j in range(0, n): - matrix[i][j] = self.basis(X1[i] - X2[j], self.parameters.radius) + matrix = cdist( + X1, X2, lambda x, y: self.basis(x-y, self.parameters.radius) + ) return matrix def _get_weights(self, X, Y): @@ -246,13 +246,12 @@ def _get_weights(self, X, Y): """ n_points = X.shape[0] dim = X.shape[1] - identity = np.ones(n_points).reshape(n_points, 1) + identity = np.ones((n_points, 1)) dist = self._distance_matrix(X, X) H = np.bmat([[dist, identity, X], [identity.T, np.zeros((1, 1)), np.zeros((1, dim))], \ [X.T, np.zeros((dim, 1)), np.zeros((dim, dim))]]) rhs = np.bmat([[Y], [np.zeros((1, dim))], [np.zeros((dim, dim))]]) - inv_matrix_h = np.linalg.inv(H) - weights = np.dot(inv_matrix_h, rhs) + weights = np.linalg.solve(H, rhs) return weights def perform(self): @@ -264,6 +263,6 @@ def perform(self): dist = self._distance_matrix( self.original_mesh_points, self.parameters.original_control_points ) - identity = np.ones(n_points).reshape(n_points, 1) + identity = np.ones((n_points, 1)) H = np.bmat([[dist, identity, self.original_mesh_points]]) self.modified_mesh_points = np.asarray(np.dot(H, self.weights)) diff --git a/tests/test_datasets/parameters_idw_default.prm b/tests/test_datasets/parameters_idw_default.prm new file mode 100644 index 0000000..5e4a5b1 --- /dev/null +++ b/tests/test_datasets/parameters_idw_default.prm @@ -0,0 +1,34 @@ + +[Inverse Distance Weighting] +# This section describes the settings of idw. + +# the power parameter +power = 2 + + +[Control points] +# This section describes the IDW control points. + +# original control points collects the coordinates +# of the interpolation control points before the +# deformation. +original control points: 0.0 0.0 0.0 + 0.0 0.0 1.0 + 0.0 1.0 0.0 + 1.0 0.0 0.0 + 0.0 1.0 1.0 + 1.0 0.0 1.0 + 1.0 1.0 0.0 + 1.0 1.0 1.0 + +# deformed control points collects the coordinates +# of the interpolation control points after the +# deformation. +deformed control points: 0.0 0.0 0.0 + 0.0 0.0 1.0 + 0.0 1.0 0.0 + 1.0 0.0 0.0 + 0.0 1.0 1.0 + 1.0 0.0 1.0 + 1.0 1.0 0.0 + 1.0 1.0 1.0 diff --git a/tests/test_datasets/parameters_idw_deform.prm b/tests/test_datasets/parameters_idw_deform.prm new file mode 100644 index 0000000..40e3843 --- /dev/null +++ b/tests/test_datasets/parameters_idw_deform.prm @@ -0,0 +1,34 @@ + +[Inverse Distance Weighting] +# This section describes the settings of idw. + +# the power parameter +power = 3 + + +[Control points] +# This section describes the IDW control points. + +# original control points collects the coordinates +# of the interpolation control points before the +# deformation. +original control points: 0.0 0.0 0.0 + 0.0 0.0 1.0 + 0.0 1.0 0.0 + 1.0 0.0 0.0 + 0.0 1.0 1.0 + 1.0 0.0 1.0 + 1.0 1.0 0.0 + 1.0 1.0 1.0 + +# deformed control points collects the coordinates +# of the interpolation control points after the +# deformation. +deformed control points: 0.0 0.0 0.0 + 0.0 0.0 1.0 + 0.0 1.0 0.0 + 1.0 0.0 0.0 + 0.0 1.0 1.0 + 1.0 0.0 1.0 + 1.0 1.0 0.0 + 1.5 1.6 1.7 diff --git a/tests/test_idw.py b/tests/test_idw.py new file mode 100644 index 0000000..3156d0a --- /dev/null +++ b/tests/test_idw.py @@ -0,0 +1,38 @@ +from unittest import TestCase +import unittest +from pygem.idw import IDW +from pygem.params_idw import IDWParameters +import numpy as np + +class TestIDW(TestCase): + def get_cube_mesh_points(self): + # Points that define a cube + nx, ny, nz = (20, 20, 20) + mesh = np.zeros((nx * ny * nz, 3)) + xv = np.linspace(0, 1, nx) + yv = np.linspace(0, 1, ny) + zv = np.linspace(0, 1, nz) + z, y, x = np.meshgrid(zv, yv, xv) + mesh = np.array([x.ravel(), y.ravel(), z.ravel()]) + original_mesh_points = mesh.T + return original_mesh_points + + def test_idw(self): + params = IDWParameters() + params.read_parameters('tests/test_datasets/parameters_idw_default.prm') + idw = IDW(params, self.get_cube_mesh_points()) + + def test_idw_perform(self): + params = IDWParameters() + params.read_parameters('tests/test_datasets/parameters_idw_default.prm') + IDW(params, self.get_cube_mesh_points()).perform() + + def test_idw_perform_deform(self): + params = IDWParameters() + expected_stretch = [ 1.19541593, 1.36081491, 1.42095073] + params.read_parameters('tests/test_datasets/parameters_idw_deform.prm') + idw = IDW(params, self.get_cube_mesh_points()) + idw.perform() + np.testing.assert_array_almost_equal( + idw.modified_mesh_points[-3], expected_stretch + ) diff --git a/tests/test_idwparams.py b/tests/test_idwparams.py new file mode 100644 index 0000000..5dd7f78 --- /dev/null +++ b/tests/test_idwparams.py @@ -0,0 +1,95 @@ +from unittest import TestCase +import unittest +from pygem.params_idw import IDWParameters +import numpy as np +import filecmp +import os + + +class TestIDWParameters(TestCase): + def test_class_members_default_p(self): + params = IDWParameters() + assert params.power == 2 + + def test_class_members_default_original_points(self): + params = IDWParameters() + cube_vertices = np.array([ + [0., 0., 0.], + [0., 0., 1.], + [0., 1., 0.], + [1., 0., 0.], + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], + [1., 1., 1.] + ]) + np.testing.assert_equal( + params.original_control_points, cube_vertices + ) + + def test_class_members_default_deformed_points(self): + params = IDWParameters() + cube_vertices = np.array([ + [0., 0., 0.], + [0., 0., 1.], + [0., 1., 0.], + [1., 0., 0.], + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], + [1., 1., 1.] + ]) + np.testing.assert_equal( + params.deformed_control_points, cube_vertices + ) + + def test_write_parameters_filename_default(self): + params = IDWParameters() + outfilename = 'parameters_rbf.prm' + outfilename_expected = 'tests/test_datasets/parameters_idw_default.prm' + params.write_parameters(outfilename) + self.assertTrue(filecmp.cmp(outfilename, outfilename_expected)) + os.remove(outfilename) + + def test_write_not_string(self): + params = IDWParameters() + with self.assertRaises(TypeError): + params.write_parameters(5) + + def test_read_deformed(self): + params = IDWParameters() + filename = 'tests/test_datasets/parameters_idw_deform.prm' + def_vertices = np.array([ + [0., 0., 0.], + [0., 0., 1.], + [0., 1., 0.], + [1., 0., 0.], + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], + [1.5, 1.6, 1.7] + ]) + params.read_parameters(filename) + np.testing.assert_equal( + params.deformed_control_points, def_vertices + ) + + def test_read_p(self): + params = IDWParameters() + filename = 'tests/test_datasets/parameters_idw_deform.prm' + params.read_parameters(filename) + assert params.power == 3 + + def test_read_not_string(self): + params = IDWParameters() + with self.assertRaises(TypeError): + params.read_parameters(5) + + def test_read_not_real_file(self): + params = IDWParameters() + with self.assertRaises(IOError): + params.read_parameters('not_real_file') + + def test_print(self): + params = IDWParameters() + print(params)