Merge branch 'master' of git://git.ucc.asn.au/ipdf/code
[ipdf/code.git] / src / tests / bezier.cpp
1 /**
2  * Tester
3  * Bezier Curves - Compute them and put on Bitmaps
4  */
5 #include "main.h"
6 #include <cmath>
7 #include <ctime> // for performance measurements
8 #include <unordered_map> // hashtable
9 #include <map>
10 #include "screen.h"
11
12 #include "pointstobitmap.h"
13
14 using namespace std;
15 using namespace IPDF;
16
17 /**
18  * Factorial
19  * Use dynamic programming / recursion
20  */
21 int Factorial(int n)
22 {
23         static unordered_map<int, int> dp;
24         static bool init = false;
25         if (!init)
26         {
27                 init = true;
28                 dp[0] = 1;
29         }
30         auto it = dp.find(n);
31         if (it != dp.end())
32                 return it->second;
33         int result = n*Factorial(n-1);
34         dp[n] = result;
35         return result;
36 }
37
38 /**
39  * Binomial coefficients
40  */
41 int BinomialCoeff(int n, int k)
42 {
43         return Factorial(n) / Factorial(k) / Factorial(n-k);
44 }
45
46 /**
47  * Bernstein Basis Polynomial
48  */
49 template <class T> T Bernstein(int k, int n, const T & u)
50 {
51         return T(BinomialCoeff(n, k)) * pow(u, k) * pow(T(1.0) - u, n-k);
52 }
53
54 /**
55  * Evaluate a Bezier at a single point u
56  */
57 template <class T> pair<T, T> Bezier(const vector<pair<T, T> > & control, const T & u)
58 {
59         pair<T,T> result(0,0);
60         for (size_t k = 0; k < control.size(); ++k)
61         {
62                 T bez(Bernstein<T>(k, control.size()-1, u));
63                 result.first += control[k].first * bez;
64                 result.second += control[k].second * bez;
65         }
66         return result;
67 }
68
69 /**
70  * Form a Bezier curve from the control points using the specified number of intervals
71  */
72 template <class T> vector<pair<T, T> > BezierCurve(const vector<pair<T, T> > & control, int intervals = 200)
73 {
74         vector<pair<T, T> > result(intervals+1);
75         T dU = T(1)/intervals;
76         T u = 0;
77         for (int i = 0; i <= intervals; ++i)
78         {
79                 result[i] = Bezier<T>(control, u);
80                 u += dU;
81         }
82         return result;
83 }
84
85
86
87 /**
88  * Make a circle using Beziers
89  * Actually just makes one Bezier for a quadrant and then mirrors it for the others
90  */
91 template <class T> vector<pair<T, T> > BezierCircle(const T & scale = 1, int intervals = 50)
92 {
93         T k = T(4) * (pow(T(2), T(0.5)) - T(1)) / T(3); // k = 4/3 * (2^(1/2) - 1) because maths says so
94         vector<pair<T, T> >control(4);
95         control[0] = pair<T,T>(0,scale);
96         control[1] = pair<T,T>(k*scale, scale);
97         control[2] = pair<T,T>(scale, k*scale);
98         control[3] = pair<T,T>(scale, 0);
99
100         auto points = BezierCurve<T>(control, intervals);
101         
102         points.reserve(intervals*4);
103         for (int i = 0; i < intervals; ++i)
104         {
105                 points.push_back(pair<T,T>(-points[i].first, points[i].second));
106         }
107         for (int i = 0; i < intervals; ++i)
108         {
109                 points.push_back(pair<T,T>(points[i].first, -points[i].second));
110         }
111         for (int i = 0; i < intervals; ++i)
112         {
113                 points.push_back(pair<T,T>(-points[i].first, -points[i].second));
114         }
115         
116         return points;
117 }
118
119 /**
120  * Test Bezier Curve with a scale applied to the *control* points
121  */
122 template <class T> vector<pair<T,T> > TestCurve(const T & scale = 1, long intervals = 50)
123 {
124         vector<pair<T,T> > control(3);
125         control[0] = pair<T,T>(0,0);
126         control[1] = pair<T,T>(0.4*scale, 0.8*scale);
127         control[2] = pair<T,T>(scale,scale);
128         return BezierCurve<T>(control, intervals);
129 }
130
131
132 /**
133  * Test the Beziers over a range of scales. Can print the Nth bezier to stdout
134  */
135 void TestBeziers(long intervals = 50, long double start = 1e-38, long double end = 1e38, int print = -1)
136 {
137         int count = 0;
138         // Choose our Bezier here
139         //#define TESTBEZIER TestCurve
140         #define TESTBEZIER BezierCircle
141
142         for (auto scale = start; scale < end; scale *= 10)
143         {
144                 auto f = TESTBEZIER<float>(scale, intervals);
145                 auto d = TESTBEZIER<double>(scale, intervals);
146                 auto l = TESTBEZIER<long double>(scale, intervals);
147
148                 // Make bitmap(s)
149                 stringstream s;
150                 s << "bezier" << count << "_" << scale << ".bmp";
151                 PointsToBitmap<float>(f, scale, 5*(count++),0,0,s.str().c_str(),false);
152
153                 #if REAL > REAL_LONG_DOUBLE
154                         auto r = TESTBEZIER<Real>(scale, intervals);
155                 #endif
156                 if (count++ == print)
157                 {
158                         for (auto i = 0; i < f.size(); ++i)
159                         {
160                                 printf("%d\t%.50lf\t%.50lf\t%.50llf\t%.50llf", f[i].first, f[i].second, d[i].first, d[i].second, l[i].first, l[i].second);
161                                 #if REAL > REAL_LONG_DOUBLE
162                                         printf("\t%.50lf\t%.50lf", r[i].first, r[i].second);
163                                 #endif  
164                                 printf("\n");
165                         }
166                 }
167
168                 // optional: Overlay images instead
169         }
170 }
171
172 /**
173  * main
174  */
175 int main(int argc, char ** argv)
176 {       
177         TestBeziers(200);
178         // Makes a movie, not really that exciting though
179         //system("for i in *.bmp; do convert $i -filter Point -resize x600 +antialias $i; done");
180         //system("ffmpeg -i bezier%d.bmp -framerate 2 -vcodec mpeg4 bezier.mp4");
181 }
182

UCC git Repository :: git.ucc.asn.au