// Thread operations
void worker(double x_start, double x_end, int num_bins, double * result)
{
	// Worker thread variables
    double local_result = 0.0;
    double local_start, local_end, local_interval, x;
    int i, local_bins;
    int rank = omp_get_thread_num();
    int num_threads = omp_get_num_threads();

    local_interval = (x_end - x_start)/num_bins;
    local_bins = num_bins / num_threads; // integer division
	local_start = x_start + rank*local_bins*local_interval;
	local_end = local_start + local_bins*local_interval;
	if (rank == (num_threads-1)) // you are last
	{
		local_end = x_end;
		local_bins = local_bins + (num_bins - (local_bins*num_threads));
	}
	x = 0.0;

	for (i = 0; i < local_bins; i++)
	{
		// 0.5 * (b1 + b2) * h
		double start, end;
		start = local_start + (local_interval*i);
		end = start + local_interval;
		local_result += 0.5*( curve(start) + curve(end) ) * local_interval;
	}
	printf("Thread: %d, interval: %lf, bins: %d, start: %lf, end: %lf, result: %lf\n", \
			rank, local_interval, local_bins, local_start, local_end, local_result);

#	pragma omp critical
	*result += local_result;
}

static inline v4su sign_test(v4sf n, struct ray ray, float a)
{
	m34sf p = ray_point(n, ray);
	v4sf v = curve(p, a);
	m34sf p0 = ray_point(n + v4sf_set1(-0.0001), ray);
	m34sf p1 = ray_point(n + v4sf_set1(0.0001), ray);
	v4sf v0 = curve(p0, a);
	v4sf v1 = curve(p1, a);
	return sign_change(v, v0) | sign_change(v, v1);
}

void Spline::DecodeFromAss(std::string const& str) {
	// Clear current
	clear();
	std::vector<float> stack;

	// Prepare
	char command = 'm';
	Vector2D pt(0, 0);

	// Tokenize the string
	boost::char_separator<char> sep(" ");
	for (auto const& token : boost::tokenizer<boost::char_separator<char>>(str, sep)) {
		double n;
		if (agi::util::try_parse(token, &n)) {
			stack.push_back(n);

			// Move
			if (stack.size() == 2 && command == 'm') {
				pt = FromScript(Vector2D(stack[0], stack[1]));
				stack.clear();

				push_back(pt);
			}

			// Line
			if (stack.size() == 2 && command == 'l') {
				SplineCurve curve(pt, FromScript(Vector2D(stack[0], stack[1])));
				push_back(curve);

				pt = curve.p2;
				stack.clear();
			}

			// Bicubic
			else if (stack.size() == 6 && command == 'b') {
				SplineCurve curve(pt,
					FromScript(Vector2D(stack[0], stack[1])),
					FromScript(Vector2D(stack[2], stack[3])),
					FromScript(Vector2D(stack[4], stack[5])));
				push_back(curve);

				pt = curve.p4;
				stack.clear();
			}
		}
		// Got something else
		else if (token.size() == 1) {
			command = token[0];
			stack.clear();
		}
	}

	if (!empty() && front().type != SplineCurve::POINT)
		push_front(pt);
}

void tst_QEasingCurve::type()
{
    {
    QEasingCurve curve(QEasingCurve::Linear);
    QCOMPARE(curve.period(), 0.3);
    QCOMPARE(curve.amplitude(), 1.0);

    curve.setPeriod(5);
    curve.setAmplitude(3);
    QCOMPARE(curve.period(), 5.0);
    QCOMPARE(curve.amplitude(), 3.0);

    curve.setType(QEasingCurve::InElastic);
    QCOMPARE(curve.period(), 5.0);
    QCOMPARE(curve.amplitude(), 3.0);
    }

    {
    QEasingCurve curve(QEasingCurve::InElastic);
    QCOMPARE(curve.period(), 0.3);
    QCOMPARE(curve.amplitude(), 1.0);
    curve.setAmplitude(2);
    QCOMPARE(curve.type(), QEasingCurve::InElastic);
    curve.setType(QEasingCurve::Linear);
    }

    {
    // check bounaries
    QEasingCurve curve(QEasingCurve::InCubic);
    QTest::ignoreMessage(QtWarningMsg, "QEasingCurve: Invalid curve type 9999");
    curve.setType((QEasingCurve::Type)9999);
    QCOMPARE(curve.type(), QEasingCurve::InCubic);
    QTest::ignoreMessage(QtWarningMsg, "QEasingCurve: Invalid curve type -9999");
    curve.setType((QEasingCurve::Type)-9999);
    QCOMPARE(curve.type(), QEasingCurve::InCubic);
    QTest::ignoreMessage(QtWarningMsg, QString::fromAscii("QEasingCurve: Invalid curve type %1")
                        .arg(QEasingCurve::NCurveTypes).toLatin1().constData());
    curve.setType(QEasingCurve::NCurveTypes);
    QCOMPARE(curve.type(), QEasingCurve::InCubic);
    QTest::ignoreMessage(QtWarningMsg, QString::fromAscii("QEasingCurve: Invalid curve type %1")
                        .arg(QEasingCurve::Custom).toLatin1().constData());
    curve.setType(QEasingCurve::Custom);
    QCOMPARE(curve.type(), QEasingCurve::InCubic);
    QTest::ignoreMessage(QtWarningMsg, QString::fromAscii("QEasingCurve: Invalid curve type %1")
                        .arg(-1).toLatin1().constData());
    curve.setType((QEasingCurve::Type)-1);
    QCOMPARE(curve.type(), QEasingCurve::InCubic);
    curve.setType(QEasingCurve::Linear);
    QCOMPARE(curve.type(), QEasingCurve::Linear);
    curve.setType(QEasingCurve::CosineCurve);
    QCOMPARE(curve.type(), QEasingCurve::CosineCurve);
    }
}

static inline v4sf bisect(i4sf l, struct ray ray, float a)
{
	i4sf k = l;
	m34sf p0 = ray_point(k.min, ray);
	v4sf v0 = curve(p0, a);
	for (int i = 0; i < 20; i++) {
		v4sf x = v4sf_set1(0.5) * (k.min + k.max);
		m34sf p1 = ray_point(x, ray);
		v4su test = sign_change(v0, curve(p1, a));
		k.min = v4sf_select(test, k.min, x);
		k.max = v4sf_select(test, x, k.max);
	}
	return v4sf_set1(0.5) * (k.min + k.max);
}

int main ()
{ int n = 100, i, ic;
  double fpi = 3.1415926 / 180.0, step, x;
  float  xray[100], y1ray[100], y2ray[100];

  step = 360. / (n - 1);

  for (i = 0; i < n; i++)
  { xray[i] = (float) (i * step);
    x = xray[i] * fpi;
    y1ray[i] = (float) sin (x);
    y2ray[i] = (float) cos (x);
  }

  metafl ("cons");
  scrmod ("revers");
  disini ();
  pagera ();
  complx ();
  axspos (450, 1800);
  axslen (2200, 1200);

  name   ("X-axis", "x");
  name   ("Y-axis", "y");

  labdig (-1, "x");
  ticks  (9, "x");
  ticks  (10, "y");

  titlin ("Demonstration of CURVE", 1);
  titlin ("SIN(X), COS(X)", 3);

  ic =   intrgb (0.95,0.95,0.95);
  axsbgd (ic);

  graf   (0.0, 360.0, 0.0, 90.0, -1.0, 1.0, -1.0, 0.5);
  setrgb (0.7, 0.7, 0.7);
  grid   (1, 1);

  color  ("fore");
  height (50);
  title  ();

  color  ("red");
  curve  (xray, y1ray, n);
  color  ("green");
  curve  (xray, y2ray, n);
  disfin ();
  return 0;
}

static inline v4su epsilon_test(v4sf n, struct ray ray, float a)
{
	float epsilon = 0.0000000001;
	m34sf p = ray_point(n, ray);
	v4sf v = curve(p, a);
	return v4sf_lt(v4sf_abs(v), v4sf_set1(epsilon));
}

void CurveFitter::func(qreal *p, qreal *hx, int m, int n, void *data)
{
    InternalData *iData = (InternalData*)data;
    if (iData->pxy + 2 != p)
        memcpy(iData->pxy + 2, p, sizeof(qreal) * m);
    curve(SPLINE_ORDER, iData->pxy, n / 2, hx, iData->ts);
}

	void main(void *arg)
	{
		// the shading space of physical sky should be in world space
		vector Iw(normalize(vto_world(I)));
		vector ray_dir = Iw;
		// we should always enable ground blur for Max material editor. 
		// when refraction is enabled, and IOR is 1.0, because in that 
		// case, the ray direction will be coplanar with X-Y plane 
		// (Y is 0.0), and after some transformations between internal 
		// space and local frame, the precision will lose and result 
		// in either +epsilon or -epsilon for Y component.
		scalar blur = i_horizon_blur();
		if (ray_dir.y >= blur)
		{
			result() = physicalsky_color(ray_dir);
		}
		else if (ray_dir.y <= 0.0f)
		{
			result() = i_ground_color();
		}
		else
		{
			color sky_c(physicalsky_color(ray_dir));
			color ground_c(i_ground_color());

			scalar factor = curve(ray_dir.y / blur);
			result() = ground_c * (1.0f - factor) + sky_c * factor;
		}
		out->Ci = result();
		out->Oi = color(0.0f);
	}

bool pos_on_bezier(const Vector2D& pos, double range, const ControlPoint& p1, const ControlPoint& p2, Vector2D& pOut, double& tOut) {
	assert(p1.segment_after == SEGMENT_CURVE);
	// Find intersections with the horizontal and vertical lines through p0
	// theoretically we would need to check in all directions, but this covers enough
	BezierCurve curve(p1, p2);
	double roots[6];
	UInt count;
	count  = solve_cubic(curve.a.y, curve.b.y, curve.c.y, curve.d.y - pos.y, roots);
	count += solve_cubic(curve.a.x, curve.b.x, curve.c.x, curve.d.x - pos.x, roots + count); // append intersections
	// take the best intersection point
	double bestDistSqr = std::numeric_limits<double>::max(); //infinity
	for(UInt i = 0 ; i < count ; ++i) {
		double t = roots[i];
		if (t >= 0 && t < 1) {
			Vector2D pnt = curve.pointAt(t);
			double distSqr = (pnt - pos).lengthSqr();
			if (distSqr < bestDistSqr) {
				bestDistSqr = distSqr;
				pOut = pnt;
				tOut = t;
			}
		}
	}
	return bestDistSqr <= range * range;
}

Bounds bezier_bounds(const Vector2D& origin, const Matrix2D& m, const ControlPoint& p1, const ControlPoint& p2) {
	assert(p1.segment_after == SEGMENT_CURVE);
	// Transform the control points
	Vector2D r1 = origin + p1.pos * m;
	Vector2D r2 = origin + (p1.pos + p1.delta_after)  * m;
	Vector2D r3 = origin + (p2.pos + p2.delta_before) * m;
	Vector2D r4 = origin + p2.pos * m;
	// First of all, the corners should be in the bounding box
	Bounds bounds(r1);
	bounds.update(r4);
	// Solve the derivative of the bezier curve to find its extremes
	// It's only a quadtratic equation :)
	BezierCurve curve(r1,r2,r3,r4);
	double roots[4];
	UInt count;
	count  = solve_quadratic(3*curve.a.x, 2*curve.b.x, curve.c.x, roots);
	count += solve_quadratic(3*curve.a.y, 2*curve.b.y, curve.c.y, roots + count);
	// now check them for min/max
	for (UInt i = 0 ; i < count ; ++i) {
		double t = roots[i];
		if (t >=0 && t <= 1) {
			bounds.update(curve.pointAt(t));
		}
	}
	return bounds;
}

/* >>>>>>>>>> EXA_6  <<<<<<<<<< */
void exa_6 (void)
{ int i, ny, nx;
  static float    x[2] = {3.f, 9.f}, y[2];
  static char *ctyp[8] = {"SOLID",  "DOT",   "DASH", "CHNDSH",
                          "CHNDOT", "DASHM", "DOTL", "DASHL"};

  setpag ("da4p");
  disini ();
  setvlt ("small");
  pagera ();
  center ();
  chncrv ("both");
  hwfont ();

  name   ("X-axis", "x");
  name   ("Y-axis", "y");

  titlin ("Demonstration of Curve", 1);
  titlin ("Line Styles", 3);

  graf   (0.f, 10.f, 0.f, 2.f, 0.f, 10.f, 0.f, 2.f);
  title  ();

  for (i = 1; i <= 8; i++)
  { y[0] = 9.5f - i;
    y[1] = 9.5f - i;
    ny = nyposn (y[0]);
    nx = nxposn (1.0f);
    messag (ctyp[i-1], nx, ny - 20);
    curve  (x, y, 2);
  }

  disfin ();
}

void qshape_create_tool_t::subdivide_segment( float param, qshape::triple_t *prev, qshape::triple_t *cur, qshape::triple_t *next) const
{
    if( prev->corner() && next->corner())
    {
	cur->set_corner( true);
	cur->set_broken( true);
        Imath::V2f q( ( next->p1() * param) + ( prev->p1() * ( 1.0f - param)));

	cur->set_p0( q);
	cur->set_p1( q);
	cur->set_p2( q);
    }
    else
    {
	cur->set_corner( false);
	cur->set_broken( false);

	bezier::curve_t<Imath::V2f, 3> curve( prev->p1(), prev->p2(), next->p0(), next->p1());
	bezier::curve_t<Imath::V2f, 3> span1, span2;

	bezier::split_curve( curve, param, span1, span2);

	prev->set_p2( span1.p[1]);

	cur->set_p0( span1.p[2]);
	cur->set_p1( span2.p[0]);
	cur->set_p2( span2.p[1]);

	next->set_p0( span2.p[2]);
    }
}

void EvaluateMultiPath(Robot& robot,const MultiPath& path,Real t,Config& q,Real xtol,Real contactol,int numIKIters)
{
  RobotCSpace space(robot);;
  RobotGeodesicManifold manifold(robot);
  GeneralizedCubicBezierCurve curve(&space,&manifold);
  Real duration,param;
  int seg=path.Evaluate(t,curve,duration,param,MultiPath::InterpLinear);
  if(seg < 0) seg = 0;
  if(seg >= path.sections.size()) seg = (int)path.sections.size()-1;
  curve.Eval(param,q);

  //solve for constraints
  bool solveIK = false;
  if(!path.settings.contains("resolution"))
    solveIK = true;
  else {
    Real res = path.settings.as<Real>("resolution");
    if(res > xtol) solveIK=true;
  }
  if(solveIK) {
    vector<IKGoal> ik;
    path.GetIKProblem(ik,seg);
    if(!ik.empty()) {
      swap(q,robot.q);
      robot.UpdateFrames();

      int iters=numIKIters;
      bool res=SolveIK(robot,ik,contactol,iters,0);
      if(!res) printf("Warning, couldn't solve IK problem at sec %d, time %g\n",seg,t);
      swap(q,robot.q);
    }
  }
}

void AttentionMap::CollectHeatPoints(XnSkeletonJoint eJoint, std::vector<XnPoint3D> &heat_points)
{
	History *history;
	if (GetHistoryForJoint (eJoint, &history))
	{
		if (history->IsStationary() == false)
		{
			history->GetPointsNewerThanTime (m_last_update_time, g_heat_points);
		}

		if (history->IsNearTarget())
		{
			g_bezier_ctrl_pts.clear();
			history->GetApproachCurveControlPoints(g_bezier_ctrl_pts);

			BezierCurveGen curve(g_bezier_ctrl_pts);

			XnPoint3D pt;
			while(curve.next_point(pt))
			{
				g_heat_points.push_back(pt);
			};
		}
	}
}

boost::optional<ModelObject> ReverseTranslator::translateCurveQuadraticLinear(
    const WorkspaceObject& workspaceObject)
{
    CurveQuadraticLinear curve(m_model);

    OptionalString s;
    OptionalDouble d;

    if ((s = workspaceObject.name())) {
        curve.setName(*s);
    }

    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient1Constant))) {
        curve.setCoefficient1Constant(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient2x))) {
        curve.setCoefficient2x(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient3x_POW_2))) {
        curve.setCoefficient3xPOW2(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient4y))) {
        curve.setCoefficient4y(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient5x_TIMES_y))) {
        curve.setCoefficient5xTIMESY(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::Coefficient6x_POW_2_TIMES_y))) {
        curve.setCoefficient6xPOW2TIMESY(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MinimumValueofx))) {
        curve.setMinimumValueofx(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MaximumValueofx))) {
        curve.setMaximumValueofx(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MinimumValueofy))) {
        curve.setMinimumValueofy(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MaximumValueofy))) {
        curve.setMaximumValueofy(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MinimumCurveOutput))) {
        curve.setMinimumCurveOutput(*d);
    }
    if ((d = workspaceObject.getDouble(Curve_QuadraticLinearFields::MaximumCurveOutput))) {
        curve.setMaximumCurveOutput(*d);
    }
    if ((s = workspaceObject.getString(Curve_QuadraticLinearFields::InputUnitTypeforX,false,true))) {
        curve.setInputUnitTypeforX(*s);
    }
    if ((s = workspaceObject.getString(Curve_QuadraticLinearFields::InputUnitTypeforY,false,true))) {
        curve.setInputUnitTypeforY(*s);
    }
    if ((s = workspaceObject.getString(Curve_QuadraticLinearFields::OutputUnitType,false,true))) {
        curve.setOutputUnitType(*s);
    }

    return curve;
}

boost::optional<ModelObject> ReverseTranslator::translateCurveFunctionalPressureDrop( 
    const WorkspaceObject& workspaceObject )
{
  CurveFunctionalPressureDrop curve(m_model);

  OptionalString s;
  OptionalDouble d;
  
  if ((s = workspaceObject.name())) {
    curve.setName(*s);
  }

  if ((d = workspaceObject.getDouble(Curve_Functional_PressureDropFields::Diameter))) {
    curve.setDiameter(*d);
  }
  if ((d = workspaceObject.getDouble(Curve_Functional_PressureDropFields::MinorLossCoefficient))) {
    curve.setMinorLossCoefficient(*d);
  }
  if ((d = workspaceObject.getDouble(Curve_Functional_PressureDropFields::Length))) {
    curve.setLength(*d);
  }
  if ((d = workspaceObject.getDouble(Curve_Functional_PressureDropFields::Roughness))) {
    curve.setRoughness(*d);
  }
  if ((d = workspaceObject.getDouble(Curve_Functional_PressureDropFields::FixedFrictionFactor))) {
    curve.setFixedFrictionFactor(*d);
  }

  return curve;
}

void CurveShape::draw(){
    std::vector<Point> tmp;
    std::vector<Point> path;
    SolidPolygon fill(fillTexture);

    for (int i = 0; i < this->size()-3; i+=3) {
      tmp.push_back(this->at(i));
      tmp.push_back(this->at(i+1));
      tmp.push_back(this->at(i+2));
      tmp.push_back(this->at(i+3));
      Curve curve(tmp,outlineTexture);
      path = curve.getPathPoints();
      fill.insert(fill.end(),path.begin(),path.end());
      path.clear();
      tmp.clear();
    }

    //draw fill
    fill.draw();

    //draw outline
    for (int i=0;i<fill.size();i++) {
        outlineTexture.draw(fill.at(i).getX(), fill.at(i).getY());
    }
}

SCurvedHistogram::Points SCurvedHistogram::getBSplinePoints( const Points& _points ) const
{
    Points bSplinePoints;
    point_list list;        // see bspline.h

    // Add again the first point with a higher value in order to prevent B-Spline algorithm from removing
    // the first value.
    list.add_point(new point(static_cast<float>( _points[0].first), static_cast<float>( _points[0].second * 2) ));

    // Add all the points
    for(const auto& pt : _points )
    {
        list.add_point(new point(static_cast<float>(pt.first), static_cast<float>(pt.second) ));
    }

    // Add again the last point
    list.add_point(new point(static_cast<float>( _points.back().first),
                             static_cast<float>( _points.back().second / 2 ) ));

    // Commpute the points of the B-Spline with external code from AHO (to be integrated here later).
    cat_curve curve( list );
    curve.m_precision = static_cast<int>(_points.size() * 5);
    curve.compute();

    for(int i = 0; i < curve.m_precision; ++i)
    {
        bSplinePoints.push_back( Point( curve.m_curve_point[i].x, curve.m_curve_point[i].y ) );
    }

    return bSplinePoints;
}

/**
    Accumulate a list of TopoDS_Edge objects along with their lengths.  We will use
    this repeatedly while we're rendering this text string.  We don't want to re-aquire
    this information for every point of every character.  Cache it here.  This method
    should be called once before each rendering session for the text string.
 */
void COrientationModifier::InitializeFromSketch()
{
    m_edges.clear();
    m_total_edge_length = 0.0;

    if (GetNumChildren() > 0)
    {
        std::list<TopoDS_Shape> wires;
        if (::ConvertSketchToFaceOrWire( GetFirstChild(), wires, false))
        {
            // Aggregate a list of TopoDS_Edge objects and each of their lengths.  We can
            // use this list to skip through edges that we're not interested in.  i.e. the
            // text won't sit on top of them.

            for (std::list<TopoDS_Shape>::iterator itWire = wires.begin(); itWire != wires.end(); itWire++)
            {
                TopoDS_Wire wire(TopoDS::Wire(*itWire));

                for(BRepTools_WireExplorer expEdge(TopoDS::Wire(wire)); expEdge.More(); expEdge.Next())
                {
                    TopoDS_Edge edge(TopoDS_Edge(expEdge.Current()));

                    BRepAdaptor_Curve curve(edge);
                    double edge_length = GCPnts_AbscissaPoint::Length(curve);
                    m_edges.push_back( std::make_pair(edge,edge_length) );
                    m_total_edge_length += edge_length;
                } // End for
            } // End for
        } // End if - then
    } // End if - then
}