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

Author: ndem0

docs/source/code.rst

@@ -20,4 +20,5 @@ Code Documentation
    stephandler
    utils
    gui
-
+   idw
+   params_idw

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:
+

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:

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

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

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

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))