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 #include <algorithm>
  12
  13 using namespace std;
  14 using namespace IPDF;
  15
  16 /**
  17  * Factorial
  18  * Use dynamic programming / recursion
  19  */
  20 int Factorial(int n)
  21 {
  22         static unordered_map<int, int> dp;
  23         static bool init = false;
  24         if (!init)
  25         {
  26                 init = true;
  27                 dp[0] = 1;
  28         }
  29         auto it = dp.find(n);
  30         if (it != dp.end())
  31                 return it->second;
  32         int result = n*Factorial(n-1);
  33         dp[n] = result;
  34         return result;
  35 }
  36
  37 /**
  38  * Binomial coefficients
  39  */
  40 int BinomialCoeff(int n, int k)
  41 {
  42         return Factorial(n) / Factorial(k) / Factorial(n-k);
  43 }
  44
  45 /**
  46  * Bernstein Basis Polynomial
  47  */
  48 template <class T> T Bernstein(int k, int n, const T & u)
  49 {
  50         return T(BinomialCoeff(n, k)) * pow(u, k) * pow(T(1.0) - u, n-k);
  51 }
  52
  53 /**
  54  * Evaluate a Bezier at a single point u
  55  */
  56 template <class T> pair<T, T> Bezier(const vector<pair<T, T> > & control, const T & u)
  57 {
  58         pair<T,T> result(0,0);
  59         for (size_t k = 0; k < control.size(); ++k)
  60         {
  61                 T bez(Bernstein<T>(k, control.size()-1, u));
  62                 result.first += control[k].first * bez;
  63                 result.second += control[k].second * bez;
  64         }
  65         return result;
  66 }
  67
  68 /**
  69  * Form a Bezier curve from the control points using the specified number of intervals
  70  */
  71 template <class T> vector<pair<T, T> > BezierCurve(const vector<pair<T, T> > & control, int intervals = 200)
  72 {
  73         vector<pair<T, T> > result(intervals+1);
  74         T dU = T(1)/intervals;
  75         T u = 0;
  76         for (int i = 0; i <= intervals; ++i)
  77         {
  78                 result[i] = Bezier<T>(control, u);
  79                 u += dU;
  80         }
  81         return result;
  82 }
  83
  84 /**
  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
  89  */
  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)
  91 {
  92
  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
  97
  98         // Pixel buffer
  99         unsigned char * pixels = new unsigned char[width*height*4];
 100
 101         // Overlay indicates we should load the bitmap into the pixel buffer first
 102         if (overlay)
 103         {
 104                 SDL_Surface * bmp = SDL_LoadBMP(filename);
 105                 if (bmp == NULL)
 106                 {
 107                         overlay = false; // bmp (probably) doesn't exist -> not an error (we hope)
 108                 }
 109                 else
 110                 {
 111                         width = bmp->w;
 112                         height = bmp->h;
 113                         for (int i = 0; i < width*height*(bmp->format->BytesPerPixel); ++i)
 114                         {
 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);
 117                         }
 118                 }
 119         }
 120         if (!overlay)
 121         {
 122                 for (int i = 0; i < width*height*4; ++i)
 123                         pixels[i] = 0xFF; // White pixels
 124         }
 125
 126
 127         typedef long double LD; // The temporary typedef, an underappreciated technique...
 128
 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);};
 132
 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);
 138
 139         pair<LD,LD> min(left.first, bottom.second);
 140         pair<LD,LD> max(right.first, top.second);
 141
 142         // Alternately, just do this:
 143         /*
 144         pair<LD,LD> min(-scale, -scale);
 145         pair<LD,LD> max(scale, scale);
 146         */
 147
 148         // Map each point to a pixel position
 149         for (auto i = 0; i < points.size(); ++i)
 150         {
 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
 155                 // Set colour
 156                 pixels[index+ri] = r;
 157                 pixels[index+gi] = g;
 158                 pixels[index+bi] = b;
 159         }
 160
 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
 165         #else
 166                 0xff000000, 0x00ff0000, 0x0000ff00, 0x000000ff
 167         #endif //SDL_BYTEORDER
 168         );
 169
 170         if (surf == NULL)
 171                 Fatal("SDL_CreateRGBSurfaceFrom(pixels...) failed - %s", SDL_GetError());
 172         if (SDL_SaveBMP(surf, filename) != 0)
 173                 Fatal("SDL_SaveBMP failed - %s", SDL_GetError());
 174
 175         // Cleanup
 176         SDL_FreeSurface(surf);
 177         delete [] pixels;
 178 }
 179
 180 /**
 181  * Make a circle using Beziers
 182  * Actually just makes one Bezier for a quadrant and then mirrors it for the others
 183  */
 184 template <class T> vector<pair<T, T> > BezierCircle(const T & scale = 1, int intervals = 50)
 185 {
 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);
 192
 193         auto points = BezierCurve<T>(control, intervals);
 194
 195         points.reserve(intervals*4);
 196         for (int i = 0; i < intervals; ++i)
 197         {
 198                 points.push_back(pair<T,T>(-points[i].first, points[i].second));
 199         }
 200         for (int i = 0; i < intervals; ++i)
 201         {
 202                 points.push_back(pair<T,T>(points[i].first, -points[i].second));
 203         }
 204         for (int i = 0; i < intervals; ++i)
 205         {
 206                 points.push_back(pair<T,T>(-points[i].first, -points[i].second));
 207         }
 208
 209         return points;
 210 }
 211
 212 /**
 213  * Test Bezier Curve with a scale applied to the *control* points
 214  */
 215 template <class T> vector<pair<T,T> > TestCurve(const T & scale = 1, long intervals = 50)
 216 {
 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);
 222 }
 223
 224
 225 /**
 226  * Test the Beziers over a range of scales. Can print the Nth bezier to stdout
 227  */
 228 void TestBeziers(long intervals = 50, long double start = 1e-38, long double end = 1e38, int print = -1)
 229 {
 230         int count = 0;
 231         // Choose our Bezier here
 232         //#define TESTBEZIER TestCurve
 233         #define TESTBEZIER BezierCircle
 234
 235         for (auto scale = start; scale < end; scale *= 10)
 236         {
 237                 auto f = TESTBEZIER<float>(scale, intervals);
 238                 auto d = TESTBEZIER<double>(scale, intervals);
 239                 auto l = TESTBEZIER<long double>(scale, intervals);
 240
 241                 // Make bitmap(s)
 242                 stringstream s;
 243                 s << "bezier" << count << "_" << scale << ".bmp";
 244                 MakeBitmap<float>(f, scale, 5*(count++),0,0,s.str().c_str(),false);
 245
 246                 #if REAL > REAL_LONG_DOUBLE
 247                         auto r = TESTBEZIER<Real>(scale, intervals);
 248                 #endif
 249                 if (count++ == print)
 250                 {
 251                         for (auto i = 0; i < f.size(); ++i)
 252                         {
 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);
 256                                 #endif
 257                                 printf("\n");
 258                         }
 259                 }
 260
 261                 // optional: Overlay images instead
 262         }
 263 }
 264
 265 /**
 266  * main
 267  */
 268 int main(int argc, char ** argv)
 269 {
 270         TestBeziers(200);
 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");
 274 }
 275