// TBezierInterpolation.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#define RESOLUTION 32
#define POINTS_LEN 10
typedef struct
{
double x, y;
}Point2D;
Point2D Point2DAdd( Point2D *left, Point2D *right)
{
Point2D newPoint = {left->x + right->x, left->y + right->y};
return newPoint;
}
Point2D Point2DSubtract( Point2D *left, Point2D *right)
{
Point2D newPoint = {left->x - right->x, left->y - right->y};
return newPoint;
}
Point2D Point2DMultiply( Point2D *left, double v)
{
Point2D newPoint = {left->x * v, left->y * v};
return newPoint;
}
void Point2DNormalize(Point2D *point)
{
double l
= sqrt(point
->x
* point
->x
+ point
->y
* point
->y
); point->x /= l;
point->y /= l;
}
typedef struct
{
Point2D points[4];
}Segment;
void SegmentCalc(Segment *seg,double t, Point2D *p)
{
double t2 = t * t;
double t3 = t2 * t;
double nt = 1.0 - t;
double nt2 = nt * nt;
double nt3 = nt2 * nt;
p->x = nt3 * seg->points[0].x + 3.0 * t * nt2 * seg->points[1].x + 3.0 * t2 * nt * seg->points[2].x + t3 * seg->points[3].x;
p->y = nt3 * seg->points[0].y + 3.0 * t * nt2 * seg->points[1].y + 3.0 * t2 * nt * seg->points[2].y + t3 * seg->points[3].y;
}
bool CalculateSpline( Point2D values[], int valuesSize, Segment bezier[])
{
int n = valuesSize - 1;
if (valuesSize < 2)
return false;
Point2D tgL= {0.0,0.0};
Point2D tgR= {0.0,0.0};
Point2D cur= {0.0,0.0};
Point2D next = Point2DSubtract(&values[1] ,&values[0]);
Point2D tmpPoint= {0.0,0.0};
Point2DNormalize(&next);
double l1 = 0.0, l2= 0.0, tmp= 0.0, x= 0.0;
--n;
for (int i = 0; i < n; ++i)
{
bezier[i].points[0] = bezier[i].points[1] = values[i];
bezier[i].points[2] = bezier[i].points[3] = values[i + 1];
cur = next;
next = Point2DSubtract(&values[i + 2],&values[i + 1]);
Point2DNormalize(&next);
tgL = tgR;
tgR = Point2DAdd(&cur,&next);
Point2DNormalize(&tgR);
if (abs(values
[i
+ 1].
y - values
[i
].
y) < DBL_EPSILON
) {
l1 = l2 = 0.0;
}
else
{
tmp = values[i + 1].x - values[i].x;
l1
= abs(tgL.
x) > DBL_EPSILON
? tmp
/ (2.0 * tgL.
x) : 1.0; l2
= abs(tgR.
x) > DBL_EPSILON
? tmp
/ (2.0 * tgR.
x) : 1.0; }
if (abs(tgL.
x) > DBL_EPSILON
&& abs(tgR.
x) > DBL_EPSILON
) {
tmp = tgL.y / tgL.x - tgR.y / tgR.x;
if (abs(tmp
) > DBL_EPSILON
) {
x = (values[i + 1].y - tgR.y / tgR.x * values[i + 1].x - values[i].y + tgL.y / tgL.x * values[i].x) / tmp;
if (x > values[i].x && x < values[i + 1].x)
{
if (tgL.y > 0.0)
{
if (l1 > l2)
l1 = 0.0;
else
l2 = 0.0;
}
else
{
if (l1 < l2)
l1 = 0.0;
else
l2 = 0.0;
}
}
}
}
tmpPoint = Point2DMultiply(&tgL , l1);
bezier[i].points[1] = Point2DAdd(& bezier[i].points[1],&tmpPoint);
tmpPoint = Point2DMultiply(&tgR , l2);
bezier[i].points[2] = Point2DSubtract(& bezier[i].points[2],&tmpPoint);
}
l1
= abs(tgL.
x) > DBL_EPSILON
? (values
[n
+ 1].
x - values
[n
].
x) / (2.0 * tgL.
x) : 1.0;
bezier[n].points[0] = bezier[n].points[1] = values[n];
bezier[n].points[2] = bezier[n].points[3] = values[n + 1];
tmpPoint = Point2DMultiply(&tgR , l1);
bezier[n].points[1] = Point2DAdd(& bezier[n].points[1],&tmpPoint);
return true;
}
int main()
{
Point2D testValues[POINTS_LEN] = {{0.0, 0.0},
{20.0, 0.0},
{45.0, -47.0},
{53.0, 335.0},
{57.0, 26.0},
{62.0, 387.0},
{ 74.0, 104.0},
{ 89.0, 0.0},
{95.0, 100.0},
{ 100.0, 0.0} };
Segment spline[100];
Point2D p = {0.0,0.0};
CalculateSpline(testValues,POINTS_LEN, spline);
for(int i = 0; i < POINTS_LEN - 1;i++)
{
for (int j = 0; j < RESOLUTION; j++)
{
SegmentCalc(&spline[i],(double)j / (double)RESOLUTION, &p);
printf("%lf %lf\n", p.
x, p.
y); }
}
printf("%lf %lf\n", testValues
[POINTS_LEN
- 1].
x,testValues
[POINTS_LEN
- 1].
y);
return 0;
}
Ly8gVEJlemllckludGVycG9sYXRpb24uYwojaW5jbHVkZSA8c3RkaW8uaD4KI2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPHN0ZGxpYi5oPgojaW5jbHVkZSA8c3RkYm9vbC5oPgojaW5jbHVkZSA8ZmxvYXQuaD4KCiNkZWZpbmUgUkVTT0xVVElPTiAzMgojZGVmaW5lICBQT0lOVFNfTEVOIDEwCgp0eXBlZGVmIHN0cnVjdAp7CiAgICBkb3VibGUgeCwgeTsKfVBvaW50MkQ7CgoKUG9pbnQyRCBQb2ludDJEQWRkKCBQb2ludDJEICpsZWZ0LCAgUG9pbnQyRCAqcmlnaHQpCnsKICAgIFBvaW50MkQgbmV3UG9pbnQgPSB7bGVmdC0+eCArIHJpZ2h0LT54LCBsZWZ0LT55ICsgcmlnaHQtPnl9OwogICAgcmV0dXJuIG5ld1BvaW50Owp9ClBvaW50MkQgUG9pbnQyRFN1YnRyYWN0KCBQb2ludDJEICpsZWZ0LCAgUG9pbnQyRCAqcmlnaHQpCnsKICAgIFBvaW50MkQgbmV3UG9pbnQgPSB7bGVmdC0+eCAtIHJpZ2h0LT54LCBsZWZ0LT55IC0gcmlnaHQtPnl9OwogICAgcmV0dXJuIG5ld1BvaW50Owp9ClBvaW50MkQgUG9pbnQyRE11bHRpcGx5KCBQb2ludDJEICpsZWZ0LCBkb3VibGUgdikKewogICAgUG9pbnQyRCBuZXdQb2ludCA9IHtsZWZ0LT54ICogdiwgbGVmdC0+eSAqIHZ9OwogICAgcmV0dXJuIG5ld1BvaW50Owp9CgoKdm9pZCBQb2ludDJETm9ybWFsaXplKFBvaW50MkQgKnBvaW50KQp7CiAgICBkb3VibGUgbCA9IHNxcnQocG9pbnQtPnggKiBwb2ludC0+eCArIHBvaW50LT55ICogcG9pbnQtPnkpOwogICAgcG9pbnQtPnggLz0gbDsKICAgIHBvaW50LT55IC89IGw7Cn0KCgp0eXBlZGVmIHN0cnVjdAp7CiAgICBQb2ludDJEIHBvaW50c1s0XTsKfVNlZ21lbnQ7Cgp2b2lkIFNlZ21lbnRDYWxjKFNlZ21lbnQgKnNlZyxkb3VibGUgdCwgUG9pbnQyRCAqcCkKewogICAgZG91YmxlIHQyID0gdCAqIHQ7CiAgICBkb3VibGUgdDMgPSB0MiAqIHQ7CiAgICBkb3VibGUgbnQgPSAxLjAgLSB0OwogICAgZG91YmxlIG50MiA9IG50ICogbnQ7CiAgICBkb3VibGUgbnQzID0gbnQyICogbnQ7CiAgICBwLT54ID0gbnQzICogc2VnLT5wb2ludHNbMF0ueCArIDMuMCAqIHQgKiBudDIgKiBzZWctPnBvaW50c1sxXS54ICsgMy4wICogdDIgKiBudCAqIHNlZy0+cG9pbnRzWzJdLnggKyB0MyAqIHNlZy0+cG9pbnRzWzNdLng7CiAgICBwLT55ID0gbnQzICogc2VnLT5wb2ludHNbMF0ueSArIDMuMCAqIHQgKiBudDIgKiBzZWctPnBvaW50c1sxXS55ICsgMy4wICogdDIgKiBudCAqIHNlZy0+cG9pbnRzWzJdLnkgKyB0MyAqIHNlZy0+cG9pbnRzWzNdLnk7Cn0KCgpib29sIENhbGN1bGF0ZVNwbGluZSggUG9pbnQyRCB2YWx1ZXNbXSwgaW50IHZhbHVlc1NpemUsIFNlZ21lbnQgYmV6aWVyW10pCnsKICAgIGludCBuID0gdmFsdWVzU2l6ZSAtIDE7CgogICAgaWYgKHZhbHVlc1NpemUgPCAyKQogICAgICAgIHJldHVybiBmYWxzZTsKCiAgICBQb2ludDJEIHRnTD0gezAuMCwwLjB9OwogICAgUG9pbnQyRCB0Z1I9IHswLjAsMC4wfTsKICAgIFBvaW50MkQgY3VyPSB7MC4wLDAuMH07CiAgICBQb2ludDJEIG5leHQgPSBQb2ludDJEU3VidHJhY3QoJnZhbHVlc1sxXSAsJnZhbHVlc1swXSk7CiAgICBQb2ludDJEIHRtcFBvaW50PSB7MC4wLDAuMH07CiAgICBQb2ludDJETm9ybWFsaXplKCZuZXh0KTsKCiAgICBkb3VibGUgbDEgPSAwLjAsIGwyPSAwLjAsIHRtcD0gMC4wLCB4PSAwLjA7CgogICAgLS1uOwoKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbjsgKytpKQogICAgewogICAgICAgIGJlemllcltpXS5wb2ludHNbMF0gPSBiZXppZXJbaV0ucG9pbnRzWzFdID0gdmFsdWVzW2ldOwogICAgICAgIGJlemllcltpXS5wb2ludHNbMl0gPSBiZXppZXJbaV0ucG9pbnRzWzNdID0gdmFsdWVzW2kgKyAxXTsKCiAgICAgICAgY3VyID0gbmV4dDsKICAgICAgICBuZXh0ID0gUG9pbnQyRFN1YnRyYWN0KCZ2YWx1ZXNbaSArIDJdLCZ2YWx1ZXNbaSArIDFdKTsKICAgICAgICBQb2ludDJETm9ybWFsaXplKCZuZXh0KTsKCiAgICAgICAgdGdMID0gdGdSOwoKICAgICAgICB0Z1IgPSBQb2ludDJEQWRkKCZjdXIsJm5leHQpOwogICAgICAgIFBvaW50MkROb3JtYWxpemUoJnRnUik7CgogICAgICAgIGlmIChhYnModmFsdWVzW2kgKyAxXS55IC0gdmFsdWVzW2ldLnkpIDwgREJMX0VQU0lMT04pCiAgICAgICAgewogICAgICAgICAgICBsMSA9IGwyID0gMC4wOwogICAgICAgIH0KICAgICAgICBlbHNlCiAgICAgICAgewogICAgICAgICAgICB0bXAgPSB2YWx1ZXNbaSArIDFdLnggLSB2YWx1ZXNbaV0ueDsKICAgICAgICAgICAgbDEgPSBhYnModGdMLngpID4gREJMX0VQU0lMT04gPyB0bXAgLyAoMi4wICogdGdMLngpIDogMS4wOwogICAgICAgICAgICBsMiA9IGFicyh0Z1IueCkgPiBEQkxfRVBTSUxPTiA/IHRtcCAvICgyLjAgKiB0Z1IueCkgOiAxLjA7CiAgICAgICAgfQoKICAgICAgICBpZiAoYWJzKHRnTC54KSA+IERCTF9FUFNJTE9OICYmIGFicyh0Z1IueCkgPiBEQkxfRVBTSUxPTikKICAgICAgICB7CiAgICAgICAgICAgIHRtcCA9IHRnTC55IC8gdGdMLnggLSB0Z1IueSAvIHRnUi54OwogICAgICAgICAgICBpZiAoYWJzKHRtcCkgPiBEQkxfRVBTSUxPTikKICAgICAgICAgICAgewogICAgICAgICAgICAgICAgeCA9ICh2YWx1ZXNbaSArIDFdLnkgLSB0Z1IueSAvIHRnUi54ICogdmFsdWVzW2kgKyAxXS54IC0gdmFsdWVzW2ldLnkgKyB0Z0wueSAvIHRnTC54ICogdmFsdWVzW2ldLngpIC8gdG1wOwogICAgICAgICAgICAgICAgaWYgKHggPiB2YWx1ZXNbaV0ueCAmJiB4IDwgdmFsdWVzW2kgKyAxXS54KQogICAgICAgICAgICAgICAgewogICAgICAgICAgICAgICAgICAgIGlmICh0Z0wueSA+IDAuMCkKICAgICAgICAgICAgICAgICAgICB7CiAgICAgICAgICAgICAgICAgICAgICAgIGlmIChsMSA+IGwyKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgbDEgPSAwLjA7CiAgICAgICAgICAgICAgICAgICAgICAgIGVsc2UKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGwyID0gMC4wOwogICAgICAgICAgICAgICAgICAgIH0KICAgICAgICAgICAgICAgICAgICBlbHNlCiAgICAgICAgICAgICAgICAgICAgewogICAgICAgICAgICAgICAgICAgICAgICBpZiAobDEgPCBsMikKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGwxID0gMC4wOwogICAgICAgICAgICAgICAgICAgICAgICBlbHNlCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBsMiA9IDAuMDsKICAgICAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICB9CiAgICAgICAgdG1wUG9pbnQgPSBQb2ludDJETXVsdGlwbHkoJnRnTCAsIGwxKTsKICAgICAgICBiZXppZXJbaV0ucG9pbnRzWzFdID0gUG9pbnQyREFkZCgmIGJlemllcltpXS5wb2ludHNbMV0sJnRtcFBvaW50KTsKICAgICAgICB0bXBQb2ludCA9IFBvaW50MkRNdWx0aXBseSgmdGdSICwgbDIpOwogICAgICAgIGJlemllcltpXS5wb2ludHNbMl0gPSBQb2ludDJEU3VidHJhY3QoJiBiZXppZXJbaV0ucG9pbnRzWzJdLCZ0bXBQb2ludCk7CiAgICB9CgogICAgbDEgPSBhYnModGdMLngpID4gREJMX0VQU0lMT04gPyAodmFsdWVzW24gKyAxXS54IC0gdmFsdWVzW25dLngpIC8gKDIuMCAqIHRnTC54KSA6IDEuMDsKCiAgICBiZXppZXJbbl0ucG9pbnRzWzBdID0gYmV6aWVyW25dLnBvaW50c1sxXSA9IHZhbHVlc1tuXTsKICAgIGJlemllcltuXS5wb2ludHNbMl0gPSBiZXppZXJbbl0ucG9pbnRzWzNdID0gdmFsdWVzW24gKyAxXTsKICAgIHRtcFBvaW50ID0gUG9pbnQyRE11bHRpcGx5KCZ0Z1IgLCBsMSk7CiAgICBiZXppZXJbbl0ucG9pbnRzWzFdID0gUG9pbnQyREFkZCgmIGJlemllcltuXS5wb2ludHNbMV0sJnRtcFBvaW50KTsKCiAgICByZXR1cm4gdHJ1ZTsKfQoKaW50IG1haW4oKQp7CgogICAgUG9pbnQyRCB0ZXN0VmFsdWVzW1BPSU5UU19MRU5dID0ge3swLjAsIDAuMH0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB7MjAuMCwgMC4wfSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHs0NS4wLCAtNDcuMH0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB7NTMuMCwgMzM1LjB9LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgezU3LjAsIDI2LjB9LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgezYyLjAsIDM4Ny4wfSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHsgNzQuMCwgMTA0LjB9LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICB7IDg5LjAsIDAuMH0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgezk1LjAsIDEwMC4wfSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICB7IDEwMC4wLCAwLjB9IH07CiAgICBTZWdtZW50IHNwbGluZVsxMDBdOwoKICAgIFBvaW50MkQgcCA9IHswLjAsMC4wfTsKCgogICAgQ2FsY3VsYXRlU3BsaW5lKHRlc3RWYWx1ZXMsUE9JTlRTX0xFTiwgc3BsaW5lKTsKCiAgICBmb3IoaW50IGkgPSAwOyBpIDwgUE9JTlRTX0xFTiAtIDE7aSsrKQogICAgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgUkVTT0xVVElPTjsgaisrKQogICAgICAgIHsKICAgICAgICAgICAgU2VnbWVudENhbGMoJnNwbGluZVtpXSwoZG91YmxlKWogLyAoZG91YmxlKVJFU09MVVRJT04sICZwKTsKICAgICAgICAgICAgcHJpbnRmKCIlbGYgJWxmXG4iLCBwLngsIHAueSk7CiAgICAgICAgfQogICAgfQoKICAgIHByaW50ZigiJWxmICVsZlxuIiwgdGVzdFZhbHVlc1tQT0lOVFNfTEVOIC0gMV0ueCx0ZXN0VmFsdWVzW1BPSU5UU19MRU4gLSAxXS55KTsKCiAgICByZXR1cm4gMDsKfQo=