/// <summary>
 /// Refer to : http://www.codeproject.com/KB/graphics/BezierSpline.aspx
 /// Solves a tridiagonal system for one of coordinates (x or y) of first Bezier control points.
 /// </summary>
 /// <param name="rhs">Right hand side vector.</param>
 /// <param name="x">Solution vector.</param>
 void GetFirstControlPoints(
     __in const std::vector<FLOAT>& rhs,
     __out std::vector<FLOAT>& x )
 {
     ATLASSERT(rhs.size()==x.size());
     int n = rhs.size();
     std::vector<FLOAT> tmp(n);    // Temp workspace.

     FLOAT b = 2.0f;
     x[0] = rhs[0] / b;
     for (int i = 1; i < n; i++) // Decomposition and forward substitution.
     {
         tmp[i] = 1 / b;
         b = (i < n-1 ? 4.0f : 3.5f) - tmp[i];
         x[i] = (rhs[i] - x[i-1]) / b;
     }
     for (int i = 1; i < n; i++)
     {
         x[n-i-1] -= tmp[n-i] * x[n-i]; // Back substitution.
     }
 }

 /// <summary>
 /// Refer to : http://www.codeproject.com/KB/graphics/BezierSpline.aspx
 /// Get open-ended Bezier Spline Control Points.
 /// </summary>
 /// <param name="knots">Input Knot Bezier spline points.</param>
 /// <param name="firstCtrlPt">Output First Control points array of knots.size()-1 length.</param>
 /// <param name="secondCtrlPt">Output Second Control points array of knots.size()-1 length.</param>
 void GetCurveControlPoints(
     __in const std::vector<D2D1_POINT_2F>& knots,
     __out std::vector<D2D1_POINT_2F>& firstCtrlPt,
     __out std::vector<D2D1_POINT_2F>& secondCtrlPt )
 {
     ATLASSERT( (firstCtrlPt.size()==secondCtrlPt.size())
         && (knots.size()==firstCtrlPt.size()+1) );

     int n = knots.size()-1;
     ATLASSERT(n>=1);

     if (n == 1)
     {
         // Special case: Bezier curve should be a straight line.
         // 3P1 = 2P0 + P3
         firstCtrlPt[0].x = (2 * knots[0].x + knots[1].x) / 3.0f;
         firstCtrlPt[0].y = (2 * knots[0].y + knots[1].y) / 3.0f;

         // P2 = 2P1 – P0
         secondCtrlPt[0].x = 2 * firstCtrlPt[0].x - knots[0].x;
         secondCtrlPt[0].y = 2 * firstCtrlPt[0].y - knots[0].y;
         return;
     }

     // Calculate first Bezier control points
     // Right hand side vector
     std::vector<FLOAT> rhs(n);

     // Set right hand side X values
     for (int i = 1; i < (n-1); ++i)
     {
         rhs[i] = 4 * knots[i].x + 2 * knots[i+1].x;
     }
     rhs[0] = knots[0].x + 2 * knots[1].x;
     rhs[n-1] = (8 * knots[n-1].x + knots[n].x) / 2.0f;
     // Get first control points X-values
     std::vector<FLOAT> x(n);
     GetFirstControlPoints(rhs,x);

     // Set right hand side Y values
     for (int i = 1; i < (n-1); ++i)
     {
         rhs[i] = 4 * knots[i].y + 2 * knots[i+1].y;
     }
     rhs[0] = knots[0].y + 2 * knots[1].y;
     rhs[n-1] = (8 * knots[n-1].y + knots[n].y) / 2.0f;
     // Get first control points Y-values
     std::vector<FLOAT> y(n);
     GetFirstControlPoints(rhs,y);

     // Fill output arrays.
     for (int i = 0; i < n; ++i)
     {
         // First control point
         firstCtrlPt[i] = D2D1::Point2F(x[i],y[i]);
         // Second control point
         if (i < (n-1))
         {
             secondCtrlPt[i] = D2D1::Point2F(2 * knots[i+1].x - x[i+1], 2*knots[i+1].y-y[i+1]);
         }
         else
         {
             secondCtrlPt[i] = D2D1::Point2F((knots[n].x + x[n-1])/2, (knots[n].y+y[n-1])/2);
         }
     }
 }

 HRESULT CreateBezierSpline(
     __in ID2D1Factory* pD2dFactory,
     __in const std::vector<D2D1_POINT_2F>& points,
     __out ID2D1PathGeometry** ppPathGeometry )
 {
     CHECK_PTR(pD2dFactory);
     CHECK_OUTPUT_PTR(ppPathGeometry);
     ATLASSERT(points.size()>1);

     int n = points.size();
     std::vector<D2D1_POINT_2F> firstCtrlPt(n-1);
     std::vector<D2D1_POINT_2F> secondCtrlPt(n-1);
     GetCurveControlPoints(points,firstCtrlPt,secondCtrlPt);

     HRESULT hr = pD2dFactory->CreatePathGeometry(ppPathGeometry);
     CHECKHR(hr);
     if (FAILED(hr))
         return hr;

     CComPtr<ID2D1GeometrySink> spSink;
     hr = (*ppPathGeometry)->Open(&spSink);
     CHECKHR(hr);
     if (SUCCEEDED(hr))
     {
         spSink->SetFillMode(D2D1_FILL_MODE_WINDING);
         spSink->BeginFigure(points[0],D2D1_FIGURE_BEGIN_FILLED);
         for (int i=1;i<n;i++)
             spSink->AddBezier(D2D1::BezierSegment(firstCtrlPt[i-1],secondCtrlPt[i-1],points[i]));
         spSink->EndFigure(D2D1_FIGURE_END_OPEN);
         spSink->Close();
     }
     return hr;
 }

 #pragma once
 #include "stdafx.h"
 #include <Direct2DHelper.h>
 using D2D1::Point2F;
 using D2D1::SizeU;
 using D2D1::ColorF;
 using D2D1::Matrix3x2F;
 using D2D1::BezierSegment;
 using D2D1::RectF;

 #include <vector>
 using std::vector;
 #include <algorithm>
 #include <boost/math/distributions/normal.hpp>

 class CMainWindow :
     public CWindowImpl<CMainWindow,CWindow,CSimpleWinTraits>
 {
 public:
     BEGIN_MSG_MAP(CMainWindow)
         MSG_WM_PAINT(OnPaint)
         MSG_WM_ERASEBKGND(OnEraseBkgnd)
         MSG_WM_SIZE(OnSize)
         MSG_WM_CREATE(OnCreate)
         MSG_WM_DESTROY(OnDestroy)
     END_MSG_MAP()

     int OnCreate(LPCREATESTRUCT /*lpCreateStruct*/)
     {
         CreateDeviceIndependentResource();
         CreateDeviceResource();
         CreateCurve();
         return 0;
     }

     void OnDestroy()
     {
         PostQuitMessage(0);
     }

     void OnPaint(CDCHandle)
     {
         CPaintDC dc(m_hWnd);
         Render();
     }

     BOOL OnEraseBkgnd(CDCHandle dc)
     {
         return TRUE;    // we have erased the background
     }

     void OnSize(UINT /*nType*/, CSize size)
     {
         if (m_spHwndRT)
         {
             m_spHwndRT->Resize(SizeU(size.cx,size.cy));
             CreateCurve();
         }
     }

 private:
     void Render()
     {
         if (!m_spHwndRT)
             CreateDeviceResource();

         m_spHwndRT->BeginDraw();
         m_spHwndRT->Clear(ColorF(ColorF::CornflowerBlue));

         m_spHwndRT->SetTransform(Matrix3x2F::Identity());

         D2D1_SIZE_F size = m_spHwndRT->GetSize();
         FLOAT width = size.width-50, height = size.height-50;
         D2D1_MATRIX_3X2_F reflectY = Direct2DHelper::ReflectYMatrix();
         D2D1_MATRIX_3X2_F translate = Matrix3x2F::Translation(size.width/2.0f,size.height/2.0f);
         m_spHwndRT->SetTransform(reflectY*translate);

         // draw coordinate axis
         m_spSolidBrush->SetColor(ColorF(ColorF::Red));
         m_spHwndRT->DrawLine(Point2F(-width*0.5f,0),Point2F(width*0.5f,0),m_spSolidBrush,2.0f);
         m_spSolidBrush->SetColor(ColorF(ColorF::DarkGreen));
         m_spHwndRT->DrawLine(Point2F(0,-height*0.5f),Point2F(0,height*0.5f),m_spSolidBrush,2.0f);

         // draw curve
         m_spSolidBrush->SetColor(ColorF(ColorF::Blue));
         m_spHwndRT->DrawGeometry(m_spPathGeometry,m_spSolidBrush,1.0f);

         // draw point marks
         m_spSolidBrush->SetColor(ColorF(ColorF::Red));
         for (auto p=m_Points.cbegin();p!=m_Points.cend();p++)
         {
             Direct2DHelper::DrawRectPoint(m_spHwndRT,m_spSolidBrush,(*p),5.0f);
         }

         HRESULT hr = m_spHwndRT->EndDraw();
         if (hr == D2DERR_RECREATE_TARGET)
             DiscardDeviceResource();
     }

     void CreateDeviceIndependentResource()
     {
         Direct2DHelper::CreateD2D1Factory(&m_spD2dFactory);
     }

     void CreateDeviceResource()
     {
         CRect rc;
         GetClientRect(&rc);

         CHECK_PTR(m_spD2dFactory);
         IFR(m_spD2dFactory->CreateHwndRenderTarget(
             D2D1::RenderTargetProperties(),
             D2D1::HwndRenderTargetProperties(m_hWnd,SizeU(rc.Width(),rc.Height())),
             &m_spHwndRT));
         IFR(m_spHwndRT->CreateSolidColorBrush(ColorF(ColorF::Red),&m_spSolidBrush));
     }

     void DiscardDeviceResource()
     {
         m_spSolidBrush.Release();
         m_spHwndRT.Release();
     }

     void CreateCurve()
     {
         if (!m_spHwndRT)
             return;
         if (m_spPathGeometry)
         {
             m_spPathGeometry.Release();
             m_Points.clear();
         }

         const int ptCount = 100;
         D2D1_SIZE_F size = m_spHwndRT->GetSize();
         FLOAT width = size.width-50.0f, height = size.height*0.4f;

 #define SIN_CURVE   
 #ifdef SIN_CURVE    // create sin curve
         FLOAT factor = static_cast<FLOAT>(4.0f*M_PI/width);
         FLOAT x = -width*0.5f, y = 0, dx = width/ptCount;
         for (int i=0;i<ptCount+1;i++)
         {
             y = height*sin(factor*x);
             m_Points.push_back(Point2F(x,y));
             x += dx;
         }
 #else                // create normal distribute curve
         FLOAT factor = 10.0f/width;
         FLOAT x = -width*0.5f, y = 0, dx = width/ptCount;
         boost::math::normal nd;
         for (int i=0;i<ptCount+1;i++)
         {
             y = height*static_cast<FLOAT>(boost::math::pdf(nd,factor*x));
             m_Points.push_back(Point2F(x,y));
             x += dx;
         }
 #endif // SIN_CURVE

         // create Bezier spline
         Direct2DHelper::CreateBezierSpline(m_spD2dFactory,m_Points,&m_spPathGeometry);
         CHECK_PTR(m_spPathGeometry);
     }

 private:
     CComPtr<ID2D1Factory> m_spD2dFactory;
     CComPtr<ID2D1HwndRenderTarget> m_spHwndRT;
     CComPtr<ID2D1SolidColorBrush> m_spSolidBrush;
     CComPtr<ID2D1PathGeometry> m_spPathGeometry;

     vector<D2D1_POINT_2F> m_Points;
 };