From: Sam Moore Date: Thu, 8 May 2014 05:30:39 +0000 (+0800) Subject: Tester for Beziers X-Git-Url: https://git.ucc.asn.au/?a=commitdiff_plain;h=669f23e6767998dec8862ecab697fe9b8d3cd2d3;p=ipdf%2Fcode.git Tester for Beziers After a brief* exciting time in which I thought Beziers looked strange at big enough scales, they don't. Not much progress on the Lit Review front... --- diff --git a/src/tests/bezier.cpp b/src/tests/bezier.cpp new file mode 100644 index 0000000..016fa5c --- /dev/null +++ b/src/tests/bezier.cpp @@ -0,0 +1,275 @@ +/** + * Tester + * Bezier Curves - Compute them and put on Bitmaps + */ +#include "main.h" +#include +#include // for performance measurements +#include // hashtable +#include +#include "screen.h" +#include + +using namespace std; +using namespace IPDF; + +/** + * Factorial + * Use dynamic programming / recursion + */ +int Factorial(int n) +{ + static unordered_map dp; + static bool init = false; + if (!init) + { + init = true; + dp[0] = 1; + } + auto it = dp.find(n); + if (it != dp.end()) + return it->second; + int result = n*Factorial(n-1); + dp[n] = result; + return result; +} + +/** + * Binomial coefficients + */ +int BinomialCoeff(int n, int k) +{ + return Factorial(n) / Factorial(k) / Factorial(n-k); +} + +/** + * Bernstein Basis Polynomial + */ +template T Bernstein(int k, int n, const T & u) +{ + return T(BinomialCoeff(n, k)) * pow(u, k) * pow(T(1.0) - u, n-k); +} + +/** + * Evaluate a Bezier at a single point u + */ +template pair Bezier(const vector > & control, const T & u) +{ + pair result(0,0); + for (size_t k = 0; k < control.size(); ++k) + { + T bez(Bernstein(k, control.size()-1, u)); + result.first += control[k].first * bez; + result.second += control[k].second * bez; + } + return result; +} + +/** + * Form a Bezier curve from the control points using the specified number of intervals + */ +template vector > BezierCurve(const vector > & control, int intervals = 200) +{ + vector > result(intervals+1); + T dU = T(1)/intervals; + T u = 0; + for (int i = 0; i <= intervals; ++i) + { + result[i] = Bezier(control, u); + u += dU; + } + return result; +} + +/** + * Map vector of points onto a Bitmap + * Because it was easier than working out the OpenGL stuff right now + * Points will be mapped according to the bounding rectangle. + * Should probably deal with possible aspect ratio difference... or just make sure to always use a square I guess + */ +template void MakeBitmap(const vector > & 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) +{ + + int ri = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 0 : 3; + int gi = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 1 : 2; + int bi = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 2 : 1; + //int ai = (SDL_BYTEORDER == SDL_LIL_ENDIAN) ? 3 : 0; // No alpha in BMP + + // Pixel buffer + unsigned char * pixels = new unsigned char[width*height*4]; + + // Overlay indicates we should load the bitmap into the pixel buffer first + if (overlay) + { + SDL_Surface * bmp = SDL_LoadBMP(filename); + if (bmp == NULL) + { + overlay = false; // bmp (probably) doesn't exist -> not an error (we hope) + } + else + { + width = bmp->w; + height = bmp->h; + for (int i = 0; i < width*height*(bmp->format->BytesPerPixel); ++i) + { + // We're assuming the BMP was stored in the same Byteorder as we will save it + pixels[i] = *((unsigned char*)(bmp->pixels)+i); + } + } + } + if (!overlay) + { + for (int i = 0; i < width*height*4; ++i) + pixels[i] = 0xFF; // White pixels + } + + + typedef long double LD; // The temporary typedef, an underappreciated technique... + + // Named lambdas... this is getting worrying... + auto lessx = [](const pair & a, const pair & b){return (a.first < b.first);}; + auto lessy = [](const pair & a, const pair & b){return (a.second < b.second);}; + + // So I don't have to type as much here + pair left = *min_element(points.begin(), points.end(), lessx); + pair right = *max_element(points.begin(), points.end(), lessx); + pair bottom = *min_element(points.begin(), points.end(), lessy); + pair top = *max_element(points.begin(), points.end(),lessy); + + pair min(left.first, bottom.second); + pair max(right.first, top.second); + + // Alternately, just do this: + /* + pair min(-scale, -scale); + pair max(scale, scale); + */ + + // Map each point to a pixel position + for (auto i = 0; i < points.size(); ++i) + { + // Do maths with long double; this way any artefacts are more likely due to the creation of the bezier itself than this mapping operation + uint64_t x = llround(((LD)(points[i].first) - (LD)(min.first)) * (LD)(width-1)/((LD)max.first - (LD)min.first)); + uint64_t y = llround(((LD)(points[i].second) - (LD)(min.second)) * (LD)(height-1)/((LD)max.second - (LD)min.second)); + int index = 4*(y*width + x); // Get index into pixel array + // Set colour + pixels[index+ri] = r; + pixels[index+gi] = g; + pixels[index+bi] = b; + } + + // Truly a thing of beauty + SDL_Surface * surf = SDL_CreateRGBSurfaceFrom(pixels, width, height, 8*4, width*4, + #if SDL_BYTEORDER == SDL_LIL_ENDIAN + 0x000000ff, 0x0000ff00, 0x00ff0000, 0xff000000 + #else + 0xff000000, 0x00ff0000, 0x0000ff00, 0x000000ff + #endif //SDL_BYTEORDER + ); + + if (surf == NULL) + Fatal("SDL_CreateRGBSurfaceFrom(pixels...) failed - %s", SDL_GetError()); + if (SDL_SaveBMP(surf, filename) != 0) + Fatal("SDL_SaveBMP failed - %s", SDL_GetError()); + + // Cleanup + SDL_FreeSurface(surf); + delete [] pixels; +} + +/** + * Make a circle using Beziers + * Actually just makes one Bezier for a quadrant and then mirrors it for the others + */ +template vector > BezierCircle(const T & scale = 1, int intervals = 50) +{ + 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 + vector >control(4); + control[0] = pair(0,scale); + control[1] = pair(k*scale, scale); + control[2] = pair(scale, k*scale); + control[3] = pair(scale, 0); + + auto points = BezierCurve(control, intervals); + + points.reserve(intervals*4); + for (int i = 0; i < intervals; ++i) + { + points.push_back(pair(-points[i].first, points[i].second)); + } + for (int i = 0; i < intervals; ++i) + { + points.push_back(pair(points[i].first, -points[i].second)); + } + for (int i = 0; i < intervals; ++i) + { + points.push_back(pair(-points[i].first, -points[i].second)); + } + + return points; +} + +/** + * Test Bezier Curve with a scale applied to the *control* points + */ +template vector > TestCurve(const T & scale = 1, long intervals = 50) +{ + vector > control(3); + control[0] = pair(0,0); + control[1] = pair(0.4*scale, 0.8*scale); + control[2] = pair(scale,scale); + return BezierCurve(control, intervals); +} + + +/** + * Test the Beziers over a range of scales. Can print the Nth bezier to stdout + */ +void TestBeziers(long intervals = 50, long double start = 1e-38, long double end = 1e38, int print = -1) +{ + int count = 0; + // Choose our Bezier here + //#define TESTBEZIER TestCurve + #define TESTBEZIER BezierCircle + + for (auto scale = start; scale < end; scale *= 10) + { + auto f = TESTBEZIER(scale, intervals); + auto d = TESTBEZIER(scale, intervals); + auto l = TESTBEZIER(scale, intervals); + + // Make bitmap(s) + stringstream s; + s << "bezier" << count << "_" << scale << ".bmp"; + MakeBitmap(f, scale, 5*(count++),0,0,s.str().c_str(),false); + + #if REAL > REAL_LONG_DOUBLE + auto r = TESTBEZIER(scale, intervals); + #endif + if (count++ == print) + { + for (auto i = 0; i < f.size(); ++i) + { + 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); + #if REAL > REAL_LONG_DOUBLE + printf("\t%.50lf\t%.50lf", r[i].first, r[i].second); + #endif + printf("\n"); + } + } + + // optional: Overlay images instead + } +} + +/** + * main + */ +int main(int argc, char ** argv) +{ + TestBeziers(200); + // Makes a movie, not really that exciting though + //system("for i in *.bmp; do convert $i -filter Point -resize x600 +antialias $i; done"); + //system("ffmpeg -i bezier%d.bmp -framerate 2 -vcodec mpeg4 bezier.mp4"); +} +