B-Spline Surface Interpolation - Python 3 -

surfaceInterpolation.py

"""
	Written by Alex Jaeger
	Derived from http://paulbourke.net/geometry/spline/
"""

import random
import numpy as np

def outputOOGL(controlMatrix, outputMatrix):	
	print("LIST")
	
	# surface polygon
	print("{ = CQUAD")
	for i in range(outputMatrix.shape[0]-1):
		for j in range(outputMatrix.shape[1]-1):
			fmtString = "{} {} {} 1 1 1 1"
			print(fmtString.format(*outputMatrix[i, j]))
			print(fmtString.format(*outputMatrix[i, j+1]))
			print(fmtString.format(*outputMatrix[i+1, j+1]))
			print(fmtString.format(*outputMatrix[i+1, j]))
	print("}")
	
	# control skeleton polygon
	for i in range(controlMatrix.shape[0]-1):
		for j in range(controlMatrix.shape[1]-1):
			print("{ = SKEL 4 1  ")
			
			fmtString = "{} {} {}"
			print(fmtString.format(*controlMatrix[i, j]))
			print(fmtString.format(*controlMatrix[i, j+1]))
			print(fmtString.format(*controlMatrix[i+1, j+1]))
			print(fmtString.format(*controlMatrix[i+1, j]))
			print("5 0 1 2 3 0")
			print("}")

class BSplineSurface:			
	def __init__(self, controlMatrix):		
		self.cm = controlMatrix
		self.cmSize = controlMatrix.shape
		self.nI = self.cmSize[0]-1
		self.nJ = self.cmSize[1]-1
	
	"""
		Returns a matrix of size outputMatrixSize where the values 
		are a surface that fits the B Spline Surface created 
		by the controlMatrix
	"""
	def interpolate(self, omSize, orderI, orderJ):		
		i = 0
		j = 0
		
		om = np.zeros(shape=omSize, dtype=(float, 3))
		
		incI = (self.nI - orderI + 2) / float(omSize[0] - 1)
		incJ = (self.nJ - orderJ + 2) / float(omSize[1] - 1)
		
		knotsI = self._calcSplineKnots(self.nI, orderI)
		knotsJ = self._calcSplineKnots(self.nJ, orderJ)
		
		intervalI = 0
		for i in range(omSize[0]-1):
			intervalJ = 0
			for j in range(omSize[1]-1):
				cP = [0, 0, 0]
				
				for ki in range(self.cmSize[0]):
					for kj in range(self.cmSize[1]):
						bi = self._calcSplineBlend(ki, orderI, knotsI, intervalI)
						bj = self._calcSplineBlend(kj, orderJ, knotsJ, intervalJ)
						
						cP[0] += bi * bj * (self.cm[ki, kj][0])
						cP[1] += bi * bj * (self.cm[ki, kj][1])
						cP[2] += bi * bj * (self.cm[ki, kj][2])
				
				om[i, j] = cP
				intervalJ += incJ
			intervalI += incI
			
		intervalI = 0
		for i in range(omSize[0]-1):
			cP = [0, 0, 0]
			
			for ki in range(self.cmSize[0]):
				bi = self._calcSplineBlend(ki, orderI, knotsI, intervalI)
				
				cP[0] += self.cm[ki, self.nJ][0] * bi
				cP[1] += self.cm[ki, self.nJ][1] * bi
				cP[2] += self.cm[ki, self.nJ][2] * bi
			
			om[i,omSize[1]-1] = cP
			intervalI += incI

		om[omSize[0]-1, omSize[1]-1] = self.cm[self.nI, self.nJ]
		
		intervalJ = 0
		for j in range(omSize[1]-1):
			cP = [0, 0, 0]
			
			for kj in range(self.cmSize[1]):
				bj = self._calcSplineBlend(kj, orderJ, knotsJ, intervalJ)
				
				cP[0] += self.cm[self.nI, kj][0] * bj
				cP[1] += self.cm[self.nI, kj][1] * bj
				cP[2] += self.cm[self.nI, kj][2] * bj
			
			om[omSize[0]-1, j] = cP
			intervalJ += incJ
		
		om[omSize[0]-1, omSize[1]-1] = self.cm[self.nI, self.nJ]
		
		return om
	
	"""
		
	
	"""
	def _calcSplineKnots(self, numControlPoints, order):
		totalKnots = numControlPoints + order + 1
		
		knots = [0 for i in range(totalKnots)]
		for i in range(totalKnots):
			if i < order:
				knots[i] = 0
			elif i <= numControlPoints:
				knots[i] = i - order + 1
			elif i > numControlPoints:
				knots[i] = numControlPoints - order + 2
				
		return knots
		
	"""
	
	
	"""
	def _calcSplineBlend(self, k, t, u, v):
		value = 0
		
		if t == 1:
			if (u[k] <= v) and (v < u[k+1]):
				value = 1
			else:
				value = 0
		else:
			if (u[k+t-1] == u[k]) and (u[k+t] == u[k+1]):
				value = 0
			elif u[k+t-1] == u[k]:
				value = (u[k+t]-v) / (u[k+t]-u[k+1]) * self._calcSplineBlend(k+1,t-1,u,v)
			elif u[k+t] == u[k+1]:
				value = (v-u[k]) / (u[k+t-1]-u[k]) * self._calcSplineBlend(k,t-1,u,v)
			else:
				value = (v-u[k]) / (u[k+t-1]-u[k]) * self._calcSplineBlend(k,t-1,u,v) + \
						(u[k+t]-v) / (u[k+t]-u[k+1]) * self._calcSplineBlend(k+1,t-1,u,v)
		return value
		
outputMatrixSize = (100, 100)
inputGrid = np.zeros((3, 3), dtype=(float, 3))
polyOrderX = 3
polyOrderY = 3

######
# Testing Data - I need to create a random surface

inputGrid[0,0] = [-10, -10, 10]
inputGrid[0,2] = [-10, 10, 10]
inputGrid[2,0] = [10, -10, 10]
inputGrid[2,2] = [10, 10, 10]

inputGrid[1,1] = [0, 0, -100]

inputGrid[1,2] = [0, 10, -50]
inputGrid[1,0] = [0, -10, -50]
inputGrid[0,1] = [-10, 0, -50]
inputGrid[2,1] = [10, 0, -50]

######

bss = BSplineSurface(inputGrid)
surface = bss.interpolate(outputMatrixSize, polyOrderX, polyOrderY)
outputOOGL(inputGrid, surface)