PathStatistics.cs 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Linq;
  4. using System.Text;
  5. using System.Drawing.Drawing2D;
  6. using System.Drawing;
  7. #pragma warning disable
  8. namespace Svg
  9. {
  10. public class PathStatistics
  11. {
  12. private const double GqBreak_TwoPoint = 0.57735026918962573;
  13. private const double GqBreak_ThreePoint = 0.7745966692414834;
  14. private const double GqBreak_FourPoint_01 = 0.33998104358485631;
  15. private const double GqBreak_FourPoint_02 = 0.86113631159405257;
  16. private const double GqWeight_FourPoint_01 = 0.65214515486254621;
  17. private const double GqWeight_FourPoint_02 = 0.34785484513745385;
  18. private PathData _data;
  19. private double _totalLength;
  20. private List<ISegment> _segments = new List<ISegment>();
  21. public double TotalLength { get { return _totalLength; } }
  22. public PathStatistics(PathData data)
  23. {
  24. _data = data;
  25. int i = 1;
  26. _totalLength = 0;
  27. ISegment newSegment;
  28. while (i < _data.Points.Length)
  29. {
  30. switch (_data.Types[i])
  31. {
  32. case 1:
  33. newSegment = new LineSegment(_data.Points[i - 1], _data.Points[i]);
  34. i++;
  35. break;
  36. case 3:
  37. newSegment = new CubicBezierSegment(_data.Points[i - 1], _data.Points[i], _data.Points[i + 1], _data.Points[i + 2]);
  38. i+= 3;
  39. break;
  40. default:
  41. throw new NotSupportedException();
  42. }
  43. newSegment.StartOffset = _totalLength;
  44. _segments.Add(newSegment);
  45. _totalLength += newSegment.Length;
  46. }
  47. }
  48. public void LocationAngleAtOffset(double offset, out PointF point, out float angle)
  49. {
  50. _segments[BinarySearchForSegment(offset, 0, _segments.Count - 1)].LocationAngleAtOffset(offset, out point, out angle);
  51. }
  52. public bool OffsetOnPath(double offset)
  53. {
  54. var seg = _segments[BinarySearchForSegment(offset, 0, _segments.Count - 1)];
  55. offset -= seg.StartOffset;
  56. return (offset >= 0 && offset <= seg.Length);
  57. }
  58. private int BinarySearchForSegment(double offset, int first, int last)
  59. {
  60. if (last == first)
  61. {
  62. return first;
  63. }
  64. else if ((last - first) == 1)
  65. {
  66. return (offset >= _segments[last].StartOffset ? last : first);
  67. }
  68. else
  69. {
  70. var mid = (last + first) / 2;
  71. if (offset < _segments[mid].StartOffset)
  72. {
  73. return BinarySearchForSegment(offset, first, mid);
  74. }
  75. else
  76. {
  77. return BinarySearchForSegment(offset, mid, last);
  78. }
  79. }
  80. }
  81. private interface ISegment
  82. {
  83. double StartOffset { get; set; }
  84. double Length { get; }
  85. void LocationAngleAtOffset(double offset, out PointF point, out float rotation);
  86. }
  87. private class LineSegment : ISegment
  88. {
  89. private double _length;
  90. private double _rotation;
  91. private PointF _start;
  92. private PointF _end;
  93. public double StartOffset { get; set; }
  94. public double Length { get { return _length; } }
  95. public LineSegment(PointF start, PointF end)
  96. {
  97. _start = start;
  98. _end = end;
  99. _length = Math.Sqrt(Math.Pow(end.X - start.X, 2) + Math.Pow(end.Y - start.Y, 2));
  100. _rotation = Math.Atan2(end.Y - start.Y, end.X - start.X) * 180 / Math.PI;
  101. }
  102. public void LocationAngleAtOffset(double offset, out PointF point, out float rotation)
  103. {
  104. offset -= StartOffset;
  105. if (offset < 0 || offset > _length) throw new ArgumentOutOfRangeException();
  106. point = new PointF((float)(_start.X + (offset / _length) * (_end.X - _start.X)),
  107. (float)(_start.Y + (offset / _length) * (_end.Y - _start.Y)));
  108. rotation = (float)_rotation;
  109. }
  110. }
  111. private class CubicBezierSegment : ISegment
  112. {
  113. private PointF _p0;
  114. private PointF _p1;
  115. private PointF _p2;
  116. private PointF _p3;
  117. private double _length;
  118. private Func<double, double> _integral;
  119. private SortedList<double, double> _lengths = new SortedList<double, double>();
  120. public double StartOffset { get; set; }
  121. public double Length { get { return _length; } }
  122. public CubicBezierSegment(PointF p0, PointF p1, PointF p2, PointF p3)
  123. {
  124. _p0 = p0;
  125. _p1 = p1;
  126. _p2 = p2;
  127. _p3 = p3;
  128. _integral = (t) => CubicBezierArcLengthIntegrand(_p0, _p1, _p2, _p3, t);
  129. _length = GetLength(0, 1, 0.00000001f);
  130. _lengths.Add(0, 0);
  131. _lengths.Add(_length, 1);
  132. }
  133. private double GetLength(double left, double right, double epsilon)
  134. {
  135. var fullInt = GaussianQuadrature(_integral, left, right, 4);
  136. return Subdivide(left, right, fullInt, 0, epsilon);
  137. }
  138. private double Subdivide(double left, double right, double fullInt, double totalLength, double epsilon)
  139. {
  140. var mid = (left + right) / 2;
  141. var leftValue = GaussianQuadrature(_integral, left, mid, 4);
  142. var rightValue = GaussianQuadrature(_integral, mid, right, 4);
  143. if (Math.Abs(fullInt - (leftValue + rightValue)) > epsilon) {
  144. var leftSub = Subdivide(left, mid, leftValue, totalLength, epsilon / 2.0);
  145. totalLength += leftSub;
  146. AddElementToTable(mid, totalLength);
  147. return Subdivide(mid, right, rightValue, totalLength, epsilon / 2.0) + leftSub;
  148. }
  149. else
  150. {
  151. return leftValue + rightValue;
  152. }
  153. }
  154. private void AddElementToTable(double position, double totalLength)
  155. {
  156. _lengths.Add(totalLength, position);
  157. }
  158. public void LocationAngleAtOffset(double offset, out PointF point, out float rotation)
  159. {
  160. offset -= StartOffset;
  161. if (offset < 0 || offset > _length) throw new ArgumentOutOfRangeException();
  162. var t = BinarySearchForParam(offset, 0, _lengths.Count - 1);
  163. point = CubicBezierCurve(_p0, _p1, _p2, _p3, t);
  164. var deriv = CubicBezierDerivative(_p0, _p1, _p2, _p3, t);
  165. rotation = (float)(Math.Atan2(deriv.Y, deriv.X) * 180.0 / Math.PI);
  166. }
  167. private double BinarySearchForParam(double length, int first, int last)
  168. {
  169. if (last == first)
  170. {
  171. return _lengths.Values[last];
  172. }
  173. else if ((last - first) == 1)
  174. {
  175. return _lengths.Values[first] + (_lengths.Values[last] - _lengths.Values[first]) *
  176. (length - _lengths.Keys[first]) / (_lengths.Keys[last] - _lengths.Keys[first]);
  177. }
  178. else
  179. {
  180. var mid = (last + first) / 2;
  181. if (length < _lengths.Keys[mid])
  182. {
  183. return BinarySearchForParam(length, first, mid);
  184. }
  185. else
  186. {
  187. return BinarySearchForParam(length, mid, last);
  188. }
  189. }
  190. }
  191. /// <summary>
  192. /// Evaluates the integral of the function over the integral using the specified number of points
  193. /// </summary>
  194. /// <param name="func"></param>
  195. /// <param name="a"></param>
  196. /// <param name="b"></param>
  197. /// <param name="points"></param>
  198. public static double GaussianQuadrature(Func<double, double> func, double a, double b, int points)
  199. {
  200. switch (points)
  201. {
  202. case 1:
  203. return (b - a) * func.Invoke((a + b) / 2.0);
  204. case 2:
  205. return (b - a) / 2.0 * (func.Invoke((b - a) / 2.0 * -1 * GqBreak_TwoPoint + (a + b) / 2.0) +
  206. func.Invoke((b - a) / 2.0 * GqBreak_TwoPoint + (a + b) / 2.0));
  207. case 3:
  208. return (b - a) / 2.0 * (5.0 / 9 * func.Invoke((b - a) / 2.0 * -1 * GqBreak_ThreePoint + (a + b) / 2.0) +
  209. 8.0 / 9 * func.Invoke((a + b) / 2.0) +
  210. 5.0 / 9 * func.Invoke((b - a) / 2.0 * GqBreak_ThreePoint + (a + b) / 2.0));
  211. case 4:
  212. return (b - a) / 2.0 * (GqWeight_FourPoint_01 * func.Invoke((b - a) / 2.0 * -1 * GqBreak_FourPoint_01 + (a + b) / 2.0) +
  213. GqWeight_FourPoint_01 * func.Invoke((b - a) / 2.0 * GqBreak_FourPoint_01 + (a + b) / 2.0) +
  214. GqWeight_FourPoint_02 * func.Invoke((b - a) / 2.0 * -1 * GqBreak_FourPoint_02 + (a + b) / 2.0) +
  215. GqWeight_FourPoint_02 * func.Invoke((b - a) / 2.0 * GqBreak_FourPoint_02 + (a + b) / 2.0));
  216. }
  217. throw new NotSupportedException();
  218. }
  219. /// <remarks>http://en.wikipedia.org/wiki/B%C3%A9zier_curve</remarks>
  220. private PointF CubicBezierCurve(PointF p0, PointF p1, PointF p2, PointF p3, double t)
  221. {
  222. return new PointF((float)(Math.Pow(1 - t, 3) * p0.X + 3 * Math.Pow(1 - t, 2) * t * p1.X +
  223. 3 * (1 - t) * Math.Pow(t, 2) * p2.X + Math.Pow(t, 3) * p3.X),
  224. (float)(Math.Pow(1 - t, 3) * p0.Y + 3 * Math.Pow(1 - t, 2) * t * p1.Y +
  225. 3 * (1 - t) * Math.Pow(t, 2) * p2.Y + Math.Pow(t, 3) * p3.Y));
  226. }
  227. /// <remarks>http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-der.html</remarks>
  228. private PointF CubicBezierDerivative(PointF p0, PointF p1, PointF p2, PointF p3, double t)
  229. {
  230. return new PointF((float)(3 * Math.Pow(1 - t, 2) * (p1.X - p0.X) + 6 * (1 - t) * t * (p2.X - p1.X) + 3 * Math.Pow(t, 2) * (p3.X - p2.X)),
  231. (float)(3 * Math.Pow(1 - t, 2) * (p1.Y - p0.Y) + 6 * (1 - t) * t * (p2.Y - p1.Y) + 3 * Math.Pow(t, 2) * (p3.Y - p2.Y)));
  232. }
  233. private double CubicBezierArcLengthIntegrand(PointF p0, PointF p1, PointF p2, PointF p3, double t)
  234. {
  235. return Math.Sqrt(Math.Pow(3 * Math.Pow(1 - t, 2) * (p1.X - p0.X) + 6 * (1 - t) * t * (p2.X - p1.X) + 3 * Math.Pow(t, 2) * (p3.X - p2.X), 2) +
  236. Math.Pow(3 * Math.Pow(1 - t, 2) * (p1.Y - p0.Y) + 6 * (1 - t) * t * (p2.Y - p1.Y) + 3 * Math.Pow(t, 2) * (p3.Y - p2.Y), 2));
  237. }
  238. }
  239. }
  240. }
  241. #pragma warning restore