001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.lib.interpolation;
011
012import imagingbook.lib.image.ImageAccessor;
013
014public class SplineInterpolator extends PixelInterpolator {
015        private final double a; 
016        private final double b;
017
018        
019        public SplineInterpolator(ImageAccessor.Scalar ia) {
020                this(0.5, 0.0); // default is a Catmull-Rom spline
021        }
022        
023        public SplineInterpolator(double a, double b) {
024                super();
025                this.a = a;
026                this.b = b;
027        }
028        
029        
030        @Override
031        public float getInterpolatedValue(ImageAccessor.Scalar ia, double x0, double y0) {
032                final int u0 = (int) Math.floor(x0);    //use floor to handle negative coordinates too
033                final int v0 = (int) Math.floor(y0);
034                double q = 0;
035                for (int j = 0; j <= 3; j++) {
036                        int v = v0 + j - 1;
037                        double p = 0;
038                        for (int i = 0; i <= 3; i++) {
039                                int u = u0 + i - 1;
040                                float pixval = ia.getVal(u, v);
041                                p = p + pixval * w_cs(x0 - u);
042                        }
043                        q = q + p * w_cs(y0 - v);
044                }
045                return (float) q;
046        }       
047        
048        private double w_cs(double x) {
049                if (x < 0) 
050                        x = -x;
051                double w = 0;
052                if (x < 1) 
053                        w = (-6*a - 9*b + 12) * x*x*x + (6*a + 12*b - 18) * x*x - 2*b + 6;
054                else if (x < 2) 
055                        w = (-6*a - b) * x*x*x + (30*a + 6*b) * x*x + (-48*a - 12*b) * x + 24*a + 8*b;
056                return w/6;
057        }
058
059}