Welcome to the Java Programming Forums


The professional, friendly Java community. 21,500 members and growing!


The Java Programming Forums are a community of Java programmers from all around the World. Our members have a wide range of skills and they all have one thing in common: A passion to learn and code Java. We invite beginner Java programmers right through to Java professionals to post here and share your knowledge. Become a part of the community, help others, expand your knowledge of Java and enjoy talking with like minded people. Registration is quick and best of all free. We look forward to meeting you.


>> REGISTER NOW TO START POSTING


Members have full access to the forums. Advertisements are removed for registered users.

Results 1 to 1 of 1

Thread: Java Tip Nov 20, 2010 - Spline Interpolation

  1. #1
    Super Moderator helloworld922's Avatar
    Join Date
    Jun 2009
    Posts
    2,896
    Thanks
    23
    Thanked 619 Times in 561 Posts
    Blog Entries
    18

    Default Java Tip Nov 20, 2010 - Spline Interpolation

    Introduction

    It's been a while since I've done one of these, but life has been pretty busy Also, this isn't the most efficient algorithm, there are several places where it could use improvement.

    Anyways, here's something that I've been working on. It's not specific to the Java programming language, but it may be useful to you if you want to perform spline interpolation over a 1D data set.

    Difficulty: Medium-hard. While this tutorial doesn't involve any difficult Java concepts (just the basic Java syntax is enough), there is significant math content, in particular systems of equations and linear algebra. It's not necessary that you completely understand the algorithm in order to know how to use it, though.

    What is Spline Interpolation?

    Spline interpolation is the action of taking a set of points and fitting multiple piece-wise continuous functions (known as splines) to these points. The most common functions to fit are polynomials. While in general if you want higher precision you'll need higher order polynomials, I have found that you can get the best accuracy using cubic functions (due to floating point errors, higher order polynomials start getting worse performance). Cubic functions are the minimum you use to get a smooth curve over all your data. Anything else will only be piecewise smooth.

    Note that spline interpolation is not the same thing as regression curve fitting! For a given set of data, you are guaranteed to cross all the data points. Also, spline interpolation can only be used on functions, so you can't have multiple y values for a given x value.

    If you're really looking for speed, simple linear functions will provide good speed with decent accuracy.

    Math Section

    The basis of spline interpolation is solve a system of linear equations. The number of equations we need (and thus the number of points needed) is dictated by the number of coefficients we need to solve for.

    For example, take the cubic function:
    a x^3 + b x^2 + c x + d = y

    This has 4 coefficients a, b, c, and d, thus would require 4 points.

    So, first step is to write down our system of equations with all the points plugged in.

    a x1^3 + b x1^2 + c x1 + d = y1
    a x2^3 + b x2^2 + c x2 + d = y2
    a x3^3 + b x3^2 + c x3 + d = y3
    a x4^3 + b x4^2 + c x4 + d = y4

    This can easily be solved using any linear-algebra solver. Because I didn't want to solve this out by hand (and didn't have access to a CAS), I chose to implement a simple Gaussian elimination algorithm with partial pivoting.

    // the gaussian elimnation algorithm
    // I won't explain this here because it's not in the scope of this tip
    	private static double[] lin_solve(double[][] matrix)
    	{
    		double[] results = new double[matrix.length];
    		int[] order = new int[matrix.length];
    		for (int i = 0; i < order.length; ++i)
    		{
    			order[i] = i;
    		}
    		for (int i = 0; i < matrix.length; ++i)
    		{
    			// partial pivot
    			int maxIndex = i;
    			for (int j = i + 1; j < matrix.length; ++j)
    			{
    				if (Math.abs(matrix[maxIndex][i]) < Math.abs(matrix[j][i]))
    				{
    					maxIndex = j;
    				}
    			}
    			if (maxIndex != i)
    			{
    				// swap order
    				{
    					int temp = order[i];
    					order[i] = order[maxIndex];
    					order[maxIndex] = temp;
    				}
    				// swap matrix
    				for (int j = 0; j < matrix[0].length; ++j)
    				{
    					double temp = matrix[i][j];
    					matrix[i][j] = matrix[maxIndex][j];
    					matrix[maxIndex][j] = temp;
    				}
    			}
    			if (Math.abs(matrix[i][i]) < 1e-15)
    			{
    				throw new RuntimeException("Singularity detected");
    			}
    			for (int j = i + 1; j < matrix.length; ++j)
    			{
    				double factor = matrix[j][i] / matrix[i][i];
    				for (int k = i; k < matrix[0].length; ++k)
    				{
    					matrix[j][k] -= matrix[i][k] * factor;
    				}
    			}
    		}
    		for (int i = matrix.length - 1; i >= 0; --i)
    		{
    			// back substitute
    			results[i] = matrix[i][matrix.length];
    			for (int j = i + 1; j < matrix.length; ++j)
    			{
    				results[i] -= results[j] * matrix[i][j];
    			}
    			results[i] /= matrix[i][i];
    		}
    		double[] correctResults = new double[results.length];
    		for (int i = 0; i < order.length; ++i)
    		{
    			// switch the order around back to the original order
    			correctResults[order[i]] = results[i];
    		}
    		return results;
    	}

    once we have the coefficients, we can easily interpolate between [x1,x4] using our cubic function. However, how do we pick x1,x2,x3, and x4?

    The simplest method is to split our data into sets which each contain the minimum number of points needed to solve the linear system of equations.

    Take the example data set:
    x data
    [1,2,3,4,5,6,7,8,9,10]

    If we're using a cubic function, we would split the data as such:

    [1,2,3,4], [5,6,7,8], [9,10]

    What do we do with the last two points? there's no possible way we could create an exact cubic function for these last two points! There are two easy solutions to this problem: one is to throw these last two data points away, and the second is to make up 2 more "fake" data points. Because spline interpolation will always intersect the data points specified, we can add extra points outside the domain (similar to adding a "dummy portion" to make a function periodic for use with Fourier Transforms). While the second method is somewhat better, I chose to implement the first method (partially because I wanted to perform some extrapolation, too), though both methods are roughly the same difficulty.

    The code

    Here's the code for the spline interpolation. This method will take in x,y data lists, the x value you want to interpolate, and the order of the polynomial you want to use. In general order 3 polynomials (cubic functions) give you the best trade-off between mathematical accuracy and floating point truncation errors.

    public static double poly_interpolate(double[] dataX, double[] dataY, double x, int power)
    	{
    		int xIndex = 0;
    		while (xIndex < dataX.length - (1 + power + (dataX.length - 1) % power) && dataX[xIndex + power] < x)
    		{
    			xIndex += power;
    		}
     
    		double matrix[][] = new double[power + 1][power + 2];
    		for (int i = 0; i < power + 1; ++i)
    		{
    			for (int j = 0; j < power; ++j)
    			{
    				matrix[i][j] = Math.pow(dataX[xIndex + i], (power - j));
    			}
    			matrix[i][power] = 1;
    			matrix[i][power + 1] = dataY[xIndex + i];
    		}
    		double[] coefficients = lin_solve(matrix);
    		double answer = 0;
    		for (int i = 0; i < coefficients.length; ++i)
    		{
    			answer += coefficients[i] * Math.pow(x, (power - i));
    		}
    		return answer;
    	}

    Conclusion

    Hopefully this may be useful to someone who wants to perform interpolation on a data set. Spline interpolation is also not limited to 2-d functions. They can be expanded to 3D data sets, though the computation time will significantly increase.

    For practical applications, it's better to "solve out" the values of the coefficients using a CAS, then plugging in the actual values later. This would eliminate the need for the linear system solver, which won't always succeed, and also has significant truncation errors associated with it.

    Also, the look-up table algorithm provided here (searching linearly through an array) isn't the most efficient method. There are other methods (such as using a tree) which will provide significant performance boosts over O(n).

    Happy coding

  2. The Following 2 Users Say Thank You to helloworld922 For This Useful Post:

    JavaPF (November 26th, 2010), nx9105 (January 23rd, 2011)


Similar Threads

  1. Java Tip Jul 29, 2010 - Swing Console Component
    By helloworld922 in forum Java Swing Tutorials
    Replies: 6
    Last Post: Yesterday, 12:08 AM
  2. Java tip Aug 26, 2010 - Circular Buffers
    By helloworld922 in forum Java Programming Tutorials
    Replies: 0
    Last Post: August 26th, 2010, 09:33 PM
  3. Java Tip Aug 4, 2010 - How to use File Filters
    By helloworld922 in forum Java Programming Tutorials
    Replies: 0
    Last Post: August 4th, 2010, 05:08 PM
  4. Java Tip Jul 5, 2010 - [Eclipse IDE] Navigating through code
    By helloworld922 in forum Java JDK & IDE Tutorials
    Replies: 1
    Last Post: July 5th, 2010, 06:28 AM
  5. [SOLVED] Interpolation Search for Strings ??
    By george in forum Algorithms & Recursion
    Replies: 3
    Last Post: May 18th, 2010, 06:11 AM