3 * Bezier Curves - Compute them and put on Bitmaps
7 #include <ctime> // for performance measurements
8 #include <unordered_map> // hashtable
18 * Use dynamic programming / recursion
22 static unordered_map<int, int> dp;
23 static bool init = false;
32 int result = n*Factorial(n-1);
38 * Binomial coefficients
40 int BinomialCoeff(int n, int k)
42 return Factorial(n) / Factorial(k) / Factorial(n-k);
46 * Bernstein Basis Polynomial
48 template <class T> T Bernstein(int k, int n, const T & u)
50 return T(BinomialCoeff(n, k)) * pow(u, k) * pow(T(1.0) - u, n-k);
54 * Evaluate a Bezier at a single point u
56 template <class T> pair<T, T> Bezier(const vector<pair<T, T> > & control, const T & u)
58 pair<T,T> result(0,0);
59 for (size_t k = 0; k < control.size(); ++k)
61 T bez(Bernstein<T>(k, control.size()-1, u));
62 result.first += control[k].first * bez;
63 result.second += control[k].second * bez;
69 * Form a Bezier curve from the control points using the specified number of intervals
71 template <class T> vector<pair<T, T> > BezierCurve(const vector<pair<T, T> > & control, int intervals = 200)
73 vector<pair<T, T> > result(intervals+1);
74 T dU = T(1)/intervals;
76 for (int i = 0; i <= intervals; ++i)
78 result[i] = Bezier<T>(control, u);
85 * Map vector of points onto a Bitmap
86 * Because it was easier than working out the OpenGL stuff right now
87 * Points will be mapped according to the bounding rectangle.
88 * Should probably deal with possible aspect ratio difference... or just make sure to always use a square I guess
90 template <class T> void MakeBitmap(const vector<pair<T, T> > & points, const T & scale = 1, uint8_t r = 0, uint8_t g = 0, uint8_t b = 0, const char * filename = "bezier.bmp", bool overlay = false, uint64_t width = 400, uint64_t height = 400)
93 int ri = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 0 : 3;
94 int gi = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 1 : 2;
95 int bi = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 2 : 1;
96 //int ai = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 3 : 0; // No alpha in BMP
99 unsigned char * pixels = new unsigned char[width*height*4];
101 // Overlay indicates we should load the bitmap into the pixel buffer first
104 SDL_Surface * bmp = SDL_LoadBMP(filename);
107 overlay = false; // bmp (probably) doesn't exist -> not an error (we hope)
113 for (int i = 0; i < width*height*(bmp->format->BytesPerPixel); ++i)
115 // We're assuming the BMP was stored in the same Byteorder as we will save it
116 pixels[i] = *((unsigned char*)(bmp->pixels)+i);
122 for (int i = 0; i < width*height*4; ++i)
123 pixels[i] = 0xFF; // White pixels
127 typedef long double LD; // The temporary typedef, an underappreciated technique...
129 // Named lambdas... this is getting worrying...
130 auto lessx = [](const pair<T, T> & a, const pair<T, T> & b){return (a.first < b.first);};
131 auto lessy = [](const pair<T, T> & a, const pair<T, T> & b){return (a.second < b.second);};
133 // So I don't have to type as much here
134 pair<T,T> left = *min_element(points.begin(), points.end(), lessx);
135 pair<T,T> right = *max_element(points.begin(), points.end(), lessx);
136 pair<T,T> bottom = *min_element(points.begin(), points.end(), lessy);
137 pair<T,T> top = *max_element(points.begin(), points.end(),lessy);
139 pair<LD,LD> min(left.first, bottom.second);
140 pair<LD,LD> max(right.first, top.second);
142 // Alternately, just do this:
144 pair<LD,LD> min(-scale, -scale);
145 pair<LD,LD> max(scale, scale);
148 // Map each point to a pixel position
149 for (auto i = 0; i < points.size(); ++i)
151 // Do maths with long double; this way any artefacts are more likely due to the creation of the bezier itself than this mapping operation
152 uint64_t x = llround(((LD)(points[i].first) - (LD)(min.first)) * (LD)(width-1)/((LD)max.first - (LD)min.first));
153 uint64_t y = llround(((LD)(points[i].second) - (LD)(min.second)) * (LD)(height-1)/((LD)max.second - (LD)min.second));
154 int index = 4*(y*width + x); // Get index into pixel array
156 pixels[index+ri] = r;
157 pixels[index+gi] = g;
158 pixels[index+bi] = b;
161 // Truly a thing of beauty
162 SDL_Surface * surf = SDL_CreateRGBSurfaceFrom(pixels, width, height, 8*4, width*4,
163 #if SDL_BYTEORDER == SDL_LIL_ENDIAN
164 0x000000ff, 0x0000ff00, 0x00ff0000, 0xff000000
166 0xff000000, 0x00ff0000, 0x0000ff00, 0x000000ff
167 #endif //SDL_BYTEORDER
171 Fatal("SDL_CreateRGBSurfaceFrom(pixels...) failed - %s", SDL_GetError());
172 if (SDL_SaveBMP(surf, filename) != 0)
173 Fatal("SDL_SaveBMP failed - %s", SDL_GetError());
176 SDL_FreeSurface(surf);
181 * Make a circle using Beziers
182 * Actually just makes one Bezier for a quadrant and then mirrors it for the others
184 template <class T> vector<pair<T, T> > BezierCircle(const T & scale = 1, int intervals = 50)
186 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
187 vector<pair<T, T> >control(4);
188 control[0] = pair<T,T>(0,scale);
189 control[1] = pair<T,T>(k*scale, scale);
190 control[2] = pair<T,T>(scale, k*scale);
191 control[3] = pair<T,T>(scale, 0);
193 auto points = BezierCurve<T>(control, intervals);
195 points.reserve(intervals*4);
196 for (int i = 0; i < intervals; ++i)
198 points.push_back(pair<T,T>(-points[i].first, points[i].second));
200 for (int i = 0; i < intervals; ++i)
202 points.push_back(pair<T,T>(points[i].first, -points[i].second));
204 for (int i = 0; i < intervals; ++i)
206 points.push_back(pair<T,T>(-points[i].first, -points[i].second));
213 * Test Bezier Curve with a scale applied to the *control* points
215 template <class T> vector<pair<T,T> > TestCurve(const T & scale = 1, long intervals = 50)
217 vector<pair<T,T> > control(3);
218 control[0] = pair<T,T>(0,0);
219 control[1] = pair<T,T>(0.4*scale, 0.8*scale);
220 control[2] = pair<T,T>(scale,scale);
221 return BezierCurve<T>(control, intervals);
226 * Test the Beziers over a range of scales. Can print the Nth bezier to stdout
228 void TestBeziers(long intervals = 50, long double start = 1e-38, long double end = 1e38, int print = -1)
231 // Choose our Bezier here
232 //#define TESTBEZIER TestCurve
233 #define TESTBEZIER BezierCircle
235 for (auto scale = start; scale < end; scale *= 10)
237 auto f = TESTBEZIER<float>(scale, intervals);
238 auto d = TESTBEZIER<double>(scale, intervals);
239 auto l = TESTBEZIER<long double>(scale, intervals);
243 s << "bezier" << count << "_" << scale << ".bmp";
244 MakeBitmap<float>(f, scale, 5*(count++),0,0,s.str().c_str(),false);
246 #if REAL > REAL_LONG_DOUBLE
247 auto r = TESTBEZIER<Real>(scale, intervals);
249 if (count++ == print)
251 for (auto i = 0; i < f.size(); ++i)
253 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);
254 #if REAL > REAL_LONG_DOUBLE
255 printf("\t%.50lf\t%.50lf", r[i].first, r[i].second);
261 // optional: Overlay images instead
268 int main(int argc, char ** argv)
271 // Makes a movie, not really that exciting though
272 //system("for i in *.bmp; do convert $i -filter Point -resize x600 +antialias $i; done");
273 //system("ffmpeg -i bezier%d.bmp -framerate 2 -vcodec mpeg4 bezier.mp4");