Class: XRVG::Fitting

Inherits:
Object
  • Object
show all
Defined in:
lib/fitting.rb

Overview

bezier = Fitting.adaptative_compute( points )

Class Method Summary collapse

Class Method Details

.adaptative_compute(pointlist, maxerror = 0.0001, maxiter = 10, tlength = nil) ⇒ Object

adaptative computation with automatic splitting if error not low enough, or if convergence is not fast enough



78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
 
# File 'lib/fitting.rb', line 78

def Fitting.adaptative_compute( pointlist, maxerror=0.0001, maxiter=10, tlength=nil )
  parameters, tlengthtmp = Fitting.initparameters( pointlist )
  if parameters == [0] * pointlist.length
    return [Bezier.single(:vector, pointlist[0], V2D::O, pointlist[0], V2D::O),0.0]
  end
  tlength ||= tlengthtmp
  niter = 0
  bezier = nil
  while true
    bezier, coeffs = Fitting.iterate( pointlist, parameters )
    error = Fitting.error( bezier, pointlist, parameters ) / tlength
    parameters = Fitting.renormalize( bezier, coeffs, pointlist, parameters )

    # pointlist.length > 8 because matching with a bezier needs at least 4 points
    if (niter > maxiter and error > maxerror and pointlist.length > 8)
	pointlists = [pointlist[0..pointlist.length/2 - 1], pointlist[pointlist.length/2 - 1 ..-1]]
	beziers = []
	errors  = []
	pointlists.each do |subpointlist|
 subbezier, suberror = Fitting.adaptative_compute( subpointlist, maxerror, maxiter, tlength )
 beziers << subbezier
 errors  << suberror
	end
	bezier = beziers.sum
	error  = errors.max
	break
    elsif (error < maxerror || niter > maxiter)
	break
    end
    perror = error
    niter += 1
  end
  return [bezier, error]
end

.bezierpiece(a, b, c, d) ⇒ Object

compute control points from polynomial bezier representation

a, b, c, d are such as 

piece( t ) = at3 + bt2 + ct + d



41
42
43
44
45
46
47
 
# File 'lib/fitting.rb', line 41

def Fitting.bezierpiece( a, b, c, d )
  p0 = d
  p1 = p0 + c / 3.0
  p2 = p1 + c / 3.0 + b / 3.0
  p3 = p0 + c + b + a
  return [p0, p1, p2, p3]
end

.compute(pointlist, maxerror = 0.01, maxiter = 100) ⇒ Object

Base method

Given a pointlist, compute the closest matching cubic bezier curve

Result is in the form [p1, pc1, pc2, p2], with [p1, pc1, pc2, p2] V2D

maxerror is normalized with curve length. In case of good match (that is pointlist can be modelized by cubic bezier curve), result error will be under maxerror. If not, result may be above maxerror. In that case, computation is stopped because error no longer decrease, or because iteration is too long.



59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
# File 'lib/fitting.rb', line 59

def Fitting.compute( pointlist, maxerror=0.01, maxiter=100 )
  parameters, tlength = Fitting.initparameters( pointlist )
  perror = 1.0
  niter = 0
  while true
    bezier, coeffs = Fitting.iterate( pointlist, parameters )
    error = Fitting.error( bezier, pointlist, parameters ) / tlength
    parameters = Fitting.renormalize( bezier, coeffs, pointlist, parameters )
    if (error < maxerror || (error-perror).abs < 0.00001 || niter > maxiter )
	break
    end
    niter += 1
  end
  return bezier
end

.error(bezier, pointlist, parameters) ⇒ Object

error is max error between points in pointlist and points sampled from bezier with parameters



144
145
146
147
148
149
150
151
152
153
154
155
156
 
# File 'lib/fitting.rb', line 144

def Fitting.error( bezier, pointlist, parameters )
  maxerror = 0.0
  container = V2D[]
  [pointlist, parameters].forzip do |point, parameter|
    # Trace("point #{point.inspect} parameter #{parameter}")
    error = (point - bezier.point( parameter, container, :parameter )).r
    if error > maxerror
	maxerror = error
    end
  end
  # Trace("Fitting.error #{maxerror}")
  return maxerror
end

.initparameters(pointlist, parameters = nil) ⇒ Object

compute first parameter t estimated values from length between consecutive points



21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
# File 'lib/fitting.rb', line 21

def Fitting.initparameters( pointlist, parameters=nil ) #:nodoc:
  lengths = [0.0]
  pointlist.pairs do |p1, p2|
    lengths.push( lengths[-1] + (p1-p2).r )
  end
  tlength = lengths[-1]
  if not parameters
    if not tlength == 0.0
	parameters = lengths.map {|length| length / tlength}
    else
	parameters = [0.0] * pointlist.length
    end
  end
  return [parameters, tlength]
end

.iterate(pointlist, parameters) ⇒ Object

iterate method compute new bezier parameters from pointlist and previous bezier parameters

Algo comes from www.tinaja.com/glib/bezdist.pdf

TODO : optimized



163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
 
# File 'lib/fitting.rb', line 163

def Fitting.iterate( pointlist, parameters )
  p0 = pointlist[0]
  p1 = pointlist[-1]

  sumt0 = parameters.map{ |t| t**0.0 }.sum
  sumt1 = parameters.map{ |t| t**1.0 }.sum
  sumt2 = parameters.map{ |t| t**2.0 }.sum
  sumt3 = parameters.map{ |t| t**3.0 }.sum
  sumt4 = parameters.map{ |t| t**4.0 }.sum
  sumt5 = parameters.map{ |t| t**5.0 }.sum
  sumt6 = parameters.map{ |t| t**6.0 }.sum

  psumt1 = [pointlist, parameters].forzip.foreach(2).map {|point, t| point * (t**1.0) }.inject(V2D::O){|sum, item| sum + item}
  psumt2 = [pointlist, parameters].forzip.foreach(2).map {|point, t| point * (t**2.0) }.inject(V2D::O){|sum, item| sum + item}
  psumt3 = [pointlist, parameters].forzip.foreach(2).map {|point, t| point * (t**3.0) }.inject(V2D::O){|sum, item| sum + item}

  coeff11 = sumt6 - 2 * sumt4 + sumt2
  coeff12 = sumt5 - sumt4 - sumt3 + sumt2

  coeff21 = coeff12
  coeff22 = sumt4 - 2 * sumt3 + sumt2

  result1 = (p0 - p1) * (sumt4 - sumt2) - p0 * (sumt3 - sumt1) + psumt3 - psumt1
  result2 = (p0 - p1) * (sumt3 - sumt2) - p0 * (sumt2 - sumt1) + psumt2 - psumt1
  
  matrix = Matrix[ [coeff11, coeff12], [coeff21, coeff22] ]
  matrixinv = matrix.inverse
  ax, bx = (matrixinv * Vector[result1.x, result2.x])[0..-1]
  ay, by = (matrixinv * Vector[result1.y, result2.y])[0..-1]

  a = V2D[ax, ay]
  b = V2D[bx, by]
  d = p0
  c = p1- (a + b + p0)

  piece = Fitting.bezierpiece( a, b, c, d )
  return [Bezier.raw( *piece ), [a, b, c, d] ]
end

.renormalize(bezier, coeffs, pointlist, parameters) ⇒ Object

algo comes from www.tinaja.com/glib/bezdist.pdf



114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
 
# File 'lib/fitting.rb', line 114

def Fitting.renormalize( bezier, coeffs, pointlist, parameters )
  a3, a2, a1, a0 = coeffs
  dxdu = Proc.new {|u| 3.0*a3.x*u**2 + 2.0*a2.x*u + a1.x}
  dydu = Proc.new {|u| 3.0*a3.y*u**2 + 2.0*a2.y*u + a1.y}
  container = V2D[]
  z    = Proc.new {|u,p4| p = bezier.point( u, container, :parameter ); (p.x - p4.x) * dxdu.call( u ) + (p.y - p4.y) * dydu.call( u )}
  newparameters = []
  [pointlist, parameters].forzip do |point, parameter|
    u1 = parameter
    if parameter < 0.99
	u2 = parameter + 0.01
    else
	u2 = parameter - 0.01
    end
    z1 = z.call(u1,point)
    z2 = z.call(u2,point)
    if z1 == z2
	u2 += 0.01
	z2 = z.call(u2,point)
    end
    if z1 == z2
	u2 -= 0.01
	z2 = z.call(u2,point)
    end
    newparameters << (z2 * u1 - z1 * u2)/(z2-z1)
  end
  return newparameters
end