Curve Point Distribution by Simpson Rule
C++. This project is still going on; it needs for some improvements. Its purpose is to equally distribute desired number of points on a curve. To do that, I applied Simpson rule by sweeping small increments on x axis instead of using full Integral calculations. Then, I rotated this curve point distribution on x axis by 12 times (including start points) to get a 3D shape. Later on, I will work on those improvements to reduce errors.
Here, I took the function y = x^2 + 1 as an example by 21 points:
|
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 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 |
#include <iostream> #include <cmath> #include <vector> #include <iomanip> struct Point { Point(float _x, float _y, float _z) : x(_x), y(_y), z(_z) {} float x; float y; float z; }; float GetFunctionValue(float const x); float GetFunctionIntegralValue(float const x); float GetIntegratedFunctionValue(float const x); float GetApproximateCurveLength(int const n, float const deltaX, float const a, float const b); int main() { float const a = -1.0f; // lower bound float const b = 1.0f; // upper bound std::vector<Point> curvePointCoordinates; // Direct result from integral calculation float const curveLength = GetFunctionIntegralValue(b) - GetFunctionIntegralValue(a); // First point as a reference start (lower bound) float xAxisValue = a; float yAxisValue = GetFunctionValue(xAxisValue); float zAxisValue = 0.0f; curvePointCoordinates.push_back(Point(xAxisValue, yAxisValue, zAxisValue)); // Applying Simpson Rule to get approximate curve length int const n = 40; float deltaX = std::abs(b-a) / n; // where n is number of segments float const approximateCurveLength = GetApproximateCurveLength(n, deltaX, a, b); // Calculating unit length according to number of points on the curve int const pointNums = 21; float const unitLength = approximateCurveLength / (pointNums - 1); // Sweeping from lower bound to upper bound to get x axis value // First and last point calculations are excluded int size = pointNums - 2; for(int k = 0; k < size; ++k) { //int n = 40; float newDeltaX = 0.0001f; while(true) { float lowerBound = curvePointCoordinates[k].x; float upperBound = lowerBound + newDeltaX * (n-1); float newLength = GetApproximateCurveLength(n, newDeltaX, lowerBound, upperBound); float error = std::abs(unitLength - newLength); if(error <= 0.0001f || newLength >= unitLength) { xAxisValue = upperBound; yAxisValue = GetFunctionValue(xAxisValue); curvePointCoordinates.push_back(Point(xAxisValue, yAxisValue, zAxisValue)); break; } newDeltaX += 0.0001f; } } // Last point xAxisValue = b; yAxisValue = GetFunctionValue(xAxisValue); curvePointCoordinates.push_back(Point(xAxisValue, yAxisValue, zAxisValue)); /*// Printing X-Y coordinates for Wavefront obj file std::cout << "Curve length by integral = " << curveLength << std::endl; std::cout << "Curve length by Simpson Rule = " << approximateCurveLength << std::endl; std::cout << "Coordinates of the points distributed on the curve:" << std::endl; std::cout << std::fixed; std::cout << std::setprecision(4); for(int i = 0; i < pointNums; ++i) { std::cout << "v " << curvePointCoordinates[i].x << " " << curvePointCoordinates[i].y << " " << curvePointCoordinates[i].z << std::endl; }*/ // Rotating distributed points around x axis by 12 times including start points std::vector<Point> shapePointCoordinates; int const rotateNums = 12; float const arcAngle = (2.0f * 3.1416f) / rotateNums; for(int i = 0; i < pointNums; ++i) { for(int j = 0; j < rotateNums; ++j) { xAxisValue = curvePointCoordinates[i].x; yAxisValue = cosf(j * arcAngle) * curvePointCoordinates[i].y; zAxisValue = sinf(j * arcAngle) * curvePointCoordinates[i].y; shapePointCoordinates.push_back(Point(xAxisValue, yAxisValue, zAxisValue)); } } // Printing XYZ coordinates for Wavefront obj file std::cout << "Coordinates of the rotated curve points distributed on the XYZ space:" << std::endl; std::cout << std::fixed; std::cout << std::setprecision(4); int const totalNums = rotateNums * pointNums; for(int i = 0; i < totalNums; ++i) { std::cout << "v " << shapePointCoordinates[i].x << " " << shapePointCoordinates[i].y << " " << shapePointCoordinates[i].z << std::endl; } curvePointCoordinates.clear(); shapePointCoordinates.clear(); return 0; } //////////////////////////////////////////////////////////////////////////////// float GetFunctionValue(float const x) { // Taking function y = x^2 + 1 as an example return x*x + 1.0f; } //////////////////////////////////////////////////////////////////////////////// float GetFunctionIntegralValue(float const x) { // Integral(sqrt(1 + (dy/dx)^2)) to get the length of the function y = x^2 + 1 float const firstPart = 2.0f*x*std::sqrt(4.0f*x*x + 1.0f); float const absValue = std::abs(2.0f*x + std::sqrt(4.0f*x*x + 1.0f)); float const secondPart = std::log(absValue) / std::log(2.71828f); return (firstPart + secondPart) / 4.0f; } //////////////////////////////////////////////////////////////////////////////// float GetIntegratedFunctionValue(float const x) { // Function inside Integral(sqrt(1 + (dy/dx)^2)) return std::sqrt(1.0f + (4.0f*x*x)); } //////////////////////////////////////////////////////////////////////////////// float GetApproximateCurveLength(int const n, float const deltaX, float const a, float const b) { // Applying Simpson rule to get approximate curve length float totalValue = GetIntegratedFunctionValue(a); // First value in the chain int size = n/2; for(int i = 0; i < size; ++i) { totalValue += 4.0f * GetIntegratedFunctionValue(a + deltaX * (1.0f + 2.0f*i)); } size--; for(int j = 0; j < size; ++j) { totalValue += 2.0f * GetIntegratedFunctionValue(a + deltaX * (2.0f + 2.0f*j)); } totalValue += GetIntegratedFunctionValue(b); // Last value in the chain return (deltaX / 3.0f) * totalValue; } |
References:
[1] Simpson Rule