summaryrefslogtreecommitdiffstats
path: root/test/monotone_cubic_interpolation
blob: 0b8912e612b3ea6930299bc019e793beb931a9cb (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/usr/bin/env python3
# -*- python -*-

# Copyright © 2014, 2015, 2016, 2017  Mattias Andrée (m@maandree.se)
# 
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.


# Test of monotone cubic interpolation
# using the Fritsch–Carlson method.
# Does not overshoot, but regular
# cubic interpolation with linear
# replacement for overshots is better.


# Load matplotlib.pyplot,
# it can take some time so
# print information about it.
print('Loading matplotlib.pyplot...')
import matplotlib.pyplot as plot
print('Done loading matplotlib.pyplot')

# Modules used for input data generation
from math import *
from random import *


def main():
    # Create a page with graphs
    fig = plot.figure()
    
    # Add graphs
    add_graph(fig, 221, [i / 15 for i in range(16)])
    add_graph(fig, 222, [sin(6 * i / 15) for i in range(16)])
    add_graph(fig, 223, [(i / 15) ** 0.5 for i in range(16)])
    #add_graph(fig, 221, [random() for i in range(16)])
    #add_graph(fig, 222, [random() for i in range(16)])
    #add_graph(fig, 223, [random() for i in range(16)])
    #add_graph(fig, 224, [random() for i in range(16)])
    #add_graph(fig, 224, [(-1) ** (i // 2) for i in range(16)])
    add_graph(fig, 224, [i / 15 * (-1) ** (i // 2) for i in range(16)])
    
    # Show graphs
    plot.show()

def add_graph(fig, graph_pos, input_values):
    '''
    Add a graph
    
    @param  fig:Figure                The page to which to add the graph
    @param  graph_pos:int             Where to place the graph
    @param  input_values:list<float>  The input values for each point
    '''
    # Interpolate data
    output_values = interpolate(input_values)
    # Number of input points
    n = len(input_values)
    # Number of output points
    m = len(output_values)
    # Create graph
    graph = fig.add_subplot(graph_pos)
    # Plot interpolated data
    graph.plot([i / (m - 1) for i in range(m)], output_values, 'b-')
    # Plot input data
    graph.plot([i / (n - 1) for i in range(n)], input_values, 'ro')


def interpolate(small, tension = 0):
    '''
    Interpolate data
    
    @param   small:list<float>  The input values for each point
    @param   tension:float      A [0, 1] value of the interpolation tension
    @return  :list<float>       The values for each point in a scaled up version
    '''
    large = [None] * len(small) ** 2
    small_, large_ = len(small) - 1, len(large) - 1
    # Basis functions
    #h00 = lambda t : (1 + 2 * t) * (1 - t) ** 2
    h10 = lambda t : t * (1 - t) ** 2
    h01 = lambda t : t ** 2 * (3 - 2 * t)
    h11 = lambda t : t ** 2 * (t - 1)
    # Tension coefficent
    c_ = 1 - tension
    ## Interpolant selection
    # Compute the slopes of the secant
    # lines between successive points
    ds = [small[i + 1] - small[i] for i in range(small_)]
    # Initialize the tangents at every
    # data point as the average of the secants
    ms = [ds[0]] + [(ds[i - 1] + ds[i]) / 2 for i in range(1, small_)] + [ds[small_ - 1]]
    βlast = 0
    for i in range(small_):
        if ds[i] == 0:
            # Two successive values are equal, ms[i],
            # must be zero to preserve monotonicity,
            # no idea to do further work on them.
            ms[i], βlast = 0, -1
            continue
        # Look for local extremums
        α, β = ms[i] / ds[i], ms[i + 1] / ds[i]
        if (α < 0) or (βlast < 0):
            # Local extremum found,
            # ensure piecewise monotonicity
            ms[i], β = 0, -1
        elif α ** 2 + β ** 2 > 9:
            # Otherwise, prevent overshoot and ensure
            # monotonicity by restricting the (α, β)
            # vector to a circle of radius 3.
            τ = 3 / (α ** 2 + β ** 2) ** 0.5
            ms[i], ms[i + 1] = τ * α * ds[i], τ * β  * ds[i]
        βlast = β
    ## Interpolation 
    for i in range(len(large)):
        # Scaling
        j = i * small_ / large_
        # Floor, weight, ceiling
        j, w, k = int(j), j % 1, min(int(j) + 1, small_)
        # Points
        pj, pk = small[j], small[k]
        # Tangents
        mj, mk = c_ * ms[j], c_ * ms[k]
        # Interpolation
        large[i] = pj + h10(w) * mj + h01(w) * (pk - pj) + h11(w) * mk
    return large

# Plot interpolation
main()