Howard said:
You'll have to post *something*, at least! Are we supposed to guess what
the problem is? We don't even know what the correct or incorrect results
are. (To be just a "wee" bit facetious: does it fly around the room
singing "Mammy", instead of sorting your fruits into apples and oranges?
)
More seriously, if you could at least post the relevant parts of your
code, along with an explanation about what is expected and what it
actually seen, that would help us (and you).
-Howard
Ok... here the code first:
---
FILE 1: porcupine.cpp
---
#include <iostream>
#include <stdlib.h>
#include "math.h"
#include "point2D.h"
using namespace std;
// Helper to calculate the normal vector
Point2D norm_vector(Point2D point)
{
return Point2D(-point.y, point.x);
}
// Helper to calculate the norm
float norm(Point2D point)
{
return sqrt( pow(point.x, 2) + pow(point.y, 2) );
}
// Helper to calculate the vectorproduct
float vectorproduct(Point2D point_a, Point2D point_b)
{
return point_a.x * point_b.y - point_a.y * point_b.x;
}
// Helper to calculate the faculty
int faculty(int n)
{
int faculty = 1;
for(n; n>0; n--)
{
faculty = faculty * n;
}
return faculty;
}
// Helper to calculate the binomial
float binomial(int n, int i)
{
return faculty(n) / (faculty(i) * faculty(n-i));
}
// Helper to calculate Bernstein
float bernstein(float t, int n, int i)
{
return binomial(n,i) * pow(t,i) * pow( (1.0-t), (n-i) );
}
// Helper to calculate Bezier
Point2D bezier(float t, int n, Point2D points[])
{
Point2D result;
int i = 0;
for(i; i<=n; i++)
{
result.x += bernstein(t,n,i) * points
.x;
result.y += bernstein(t,n,i) * points.y;
}
return result;
}
// Helper to calculate the 1st derivative of a Bezier curve
Point2D bezier_first_derivative(float t, int n, Point2D points[])
{
Point2D result;
int intervallborder_s = 0;
int intervallborder_t = 1;
int i = 0;
for(i; i<=n-1; i++)
{
result.x += bernstein(t,n-1,i) * (points[i+1].x - points.x);
result.y += bernstein(t,n-1,i) * (points[i+1].y - points.y);
}
result.x = result.x * (n / (intervallborder_t - intervallborder_s));
result.y = result.y * (n / (intervallborder_t - intervallborder_s));
return result;
}
// Helper to calculate the 1st derivative of a Bezier curve
Point2D bezier_second_derivative(float t, int n, Point2D points[])
{
Point2D result;
int intervallborder_s = 0;
int intervallborder_t = 1;
int i = 0;
for(i; i<=n-2; i++)
{
result.x += bernstein(t,n-2,i) * (points[i+2].x - 2*points[i+1].x
+
points.x);
result.y += bernstein(t,n-2,i) * (points[i+2].y - 2*points[i+1].y
+
points.y);
}
result.x = result.x * ( n*(n-1.0) / pow( (double)(intervallborder_t -
intervallborder_s), 2) );
result.y = result.y * ( n*(n-1.0) / pow( (double)(intervallborder_t -
intervallborder_s), 2) );
return result;
}
// Calculate curvature
float curvature(float t, int n, Point2D points[])
{
Point2D bsd = bezier_second_derivative(t,n,points);
Point2D bfd = bezier_first_derivative(t,n,points);
float vp = vectorproduct(bsd, bfd);
float n_bfd = norm(bfd);
return vp / pow(n_bfd, 3);
}
// Calculate curvature vector
Point2D curvature_vector(float t, int n, Point2D points[])
{
Point2D bfd = bezier_first_derivative(t,n,points);
// Calculate normal vector
Point2D norm_bfd = norm_vector(bfd);
// Multiply with curvature
float c = curvature(t,n,points);
Point2D result = Point2D(norm_bfd.x * c, norm_bfd.y * c);
return result;
}
// Main
int main(int argc, char *argv[])
{
// Read the value specifying the degree
int n = atoi(argv[1]);
// Read the corresponding Bezier points
Point2D points[n];
int j = 0;
for (int i=2; i<(2*(n+1))+2; i+=2)
{
points[j] = Point2D(atof(argv), atof(argv[i+1]));
j += 1;
}
// Set a start value for t
float t = 0.0;
// Iterate and print the p_i's and the corresponding curvature values
for (int k=0; k<=20; k++)
{
Point2D pi = bezier(t, n, points);
//Point2D pi = curvature_vector(t, n, points);
//cout << curvature(t,n,points) << endl;
cout << pi.x << "\t\t" << pi.y << endl;
t = t + 0.05;
}
return 0;
}
---
File 2: point2d.h
---
class Point2D
{
public:
float x;
float y;
Point2D();
Point2D(float x, float y);
};
Point2D:oint2D()
{
this->x=0;
this->y=0;
}
Point2D:oint2D(float x, float y)
{
this->x=x;
this->y=y;
}
This is what I get:
bash-2.05b# ./porcupine 2 0 0 1 1 2 1
13073.4 9.21956e-41
11798.7 8.32063e-41
10589.4 7.4678e-41
9445.51 6.66107e-41
8366.96 5.90059e-41
7353.77 5.18607e-41
6405.95 4.51765e-41
5523.5 3.89533e-41
4706.41 3.31898e-41
3954.7 2.78886e-41
3268.34 2.30486e-41
2647.36 1.86695e-41
2091.74 1.47515e-41
1601.49 1.12945e-41
1176.6 8.29709e-42
817.085 5.76354e-42
522.934 3.68962e-42
294.15 2.07532e-42
130.733 9.23456e-43
32.6832 2.31214e-43
1.85784e-10 1.4013e-45
And this is what is expected (and what I get under windows):
