grbl-LPC-CoreXY/arc_algorithm/arc.rb

426 lines
11 KiB
Ruby

require "pp"
# A tiny ruby script to design and test the implementation of the arc-support
class Numeric
# -1 if the number < 0, 1 if number >= 0
def sign
return 0 if self == 0
(self < 0) ? -1 : 1
end
end
class CircleTest
include Math
def init
@pixels = []
@tool_position = [14,14]
40.times { @pixels << '.'*40 }
end
def plot_pixel(x,y, c)
return if x<0 || y<0 || x>39 || y > 39
@pixels[y] = @pixels[y][0..x][0..-2]+c+@pixels[y][(x+1)..-1]
end
def show
@pixels.each do |line|
puts line.gsub('.','. ').gsub('0','0 ').gsub('1','1 ').gsub('2','2 ').gsub('X','X ').gsub('o','o ')
end
end
# dP[x+1,y]: 1 + 2 x
# dP[x-1,y]: 1 - 2 x
# dP[x, y+1]: 1 + 2 y
# dP[x, y-1]: 1 - 2 y
# dP[x+1, y+1]: 2 (1 + x + y) 1+2x+1+2y
# dP[x+1, y-1]: 2 (1 + x - y) 1+2x+1-2y
# dP[x-1, y-1]: 2 (1 - x - y) 2-2x-2y
# dP[x-1, y+1]: 2 (1 - x + y) 2-2x+2x
# dP[x+a, y+b]: |dx| - 2*dx*x + |dy| + 2*dy*y
# Algorithm from the wikipedia aricle on the Midpoint circle algorithm.
def raster_circle(radius)
f = 1-radius
ddF_x = 1
ddF_y = -2*radius
x = 0
y = radius
while x<=y
if f>0
y -= 1
ddF_y += 2
f += ddF_y
end
x += 1
ddF_x += 2
f += ddF_x
plot_pixel(x+14,-y+14,'X')
end
x += 1
ddF_x += 2
f += ddF_x
while y>0
if f<0
x += 1
ddF_x += 2
f += ddF_x
end
y -= 1
ddF_y += 2
f += ddF_y
plot_pixel(x+14,-y+14,'o')
end
end
# An "ideal" arc. Computationally expensive, but always pure
def pure_arc(theta, segment, angular_direction, radius)
var = [(sin(theta)*(radius)), (cos(theta)*(radius))]
error_sum = 0.0
error_count = 0.0
max_error = 0.0
(PI*2*radius).ceil.times do
if var[0].abs<var[1].abs
major_axis = 0
minor_axis = 1
else
major_axis = 1
minor_axis = 0
end
delta = [var[1].sign*angular_direction, -var[0].sign*angular_direction]
delta[0] = angular_direction if delta[0] == 0
delta[1] = angular_direction if delta[1] == 0
var[major_axis] += delta[major_axis]
var[minor_axis] = (sqrt((radius**2-var[major_axis]**2))*var[minor_axis].sign).round
error_count += 1
error = (var[minor_axis].abs-(sqrt((radius**2)-var[major_axis]**2)).abs).abs
max_error = [max_error,error].max
error_sum += error
plot_pixel(var[0]+14, -var[1]+14, 'X')
end
puts "Average error: #{error_sum/error_count} Maximum eror: #{max_error}"
end
# A DDA-direct search circle interpolator. Optimal and impure
def arc_clean(theta, angular_travel, radius)
radius = radius
x = (sin(theta)*radius).round
y = (cos(theta)*radius).round
angular_direction = angular_travel.sign
tx = (sin(theta+angular_travel)*radius).round
ty = (cos(theta+angular_travel)*radius).round
f = (x**2 + y**2 - radius**2).round
min_x = 0
max_x = 0
i = 0
pp [x,y]
while true
if i > 0
plot_pixel(x+14, -y+14, "012"[i%3].chr)
else
plot_pixel(x+14, -y+14, "X")
end
dx = (y==0) ? -x.sign : y.sign*angular_direction
dy = (x==0) ? -y.sign : -x.sign*angular_direction
pp [[x,y],[dx,dy]]
if x.abs<y.abs
f_straight = f + 1+2*x*dx
f_diagonal = f_straight + 1+2*y*dy
x += dx
if (f_straight.abs < f_diagonal.abs)
f = f_straight
else
y += dy
f = f_diagonal
end
else
f_straight = f + 1+2*y*dy
f_diagonal = f_straight + 1+2*x*dx
y += dy
if (f_straight.abs < f_diagonal.abs)
f = f_straight
else
x += dx
f = f_diagonal
end
end
f_should_be = (x**2+y**2-radius**2)
if f.round != f_should_be.round
raise "f out of range. Is #{f}, should be #{f_should_be}"
end
min_x = [x,min_x].min
max_x = [x,max_x].max
break if (x.sign == tx.sign && y.sign == ty.sign) && (x.abs>=tx.abs) && (y.abs>=ty.abs)
i += 1
end
puts "Target #{[tx,ty].inspect}"
plot_pixel(tx+14, -ty+14, "o")
pp [x,y]
puts "Diameter: #{max_x-min_x}"
end
# A DDA-direct search circle interpolator. Optimal and impure
def arc_supaoptimal(theta, angular_travel, radius)
radius = radius
x = (sin(theta)*radius).round
y = (cos(theta)*radius).round
angular_direction = angular_travel.sign
tx = (sin(theta+angular_travel)*(radius-0.5)).floor
ty = (cos(theta+angular_travel)*(radius-0.5)).floor
f = (x**2 + y**2 - radius**2).round
x2 = 2*x
y2 = 2*y
dx = (y==0) ? -x.sign : y.sign*angular_direction
dy = (x==0) ? -y.sign : -x.sign*angular_direction
max_steps = (angular_travel.abs*radius*2).floor
# dP[x+1,y]: 1 + 2 x
# dP[x, y+1]: 1 + 2 y
max_steps.times do |i|
if i > 0
plot_pixel(x+20, -y+20, "012"[i%3].chr)
else
plot_pixel(x+20, -y+20, "X")
end
raise "x2 out of range" unless x2 == 2*x
raise "y2 out of range" unless y2 == 2*y
f_should_be = (x**2+y**2-radius**2)
if f.round != f_should_be.round
show
raise "f out of range. Is #{f}, should be #{f_should_be}"
end
if x.abs<y.abs
x += dx
f += 1+x2*dx
x2 += 2*dx
f_diagonal = f + 1 + y2*dy
if (f.abs >= f_diagonal.abs)
y += dy
dx = y.sign*angular_direction unless y == 0
y2 += 2*dy
f = f_diagonal
end
dy = -x.sign*angular_direction unless x == 0
else
y += dy
f += 1+y2*dy
y2 += 2*dy
f_diagonal = f + 1 + x2*dx
if (f.abs >= f_diagonal.abs)
x += dx
dy = -x.sign*angular_direction unless x == 0
x2 += 2*dx
f = f_diagonal
end
dx = y.sign*angular_direction unless y == 0
end
break if x*ty.sign*angular_direction>=tx*ty.sign*angular_direction && y*tx.sign*angular_direction<=ty*tx.sign*angular_direction
end
plot_pixel(tx+20, -ty+20, "o")
return {:tx => tx, :ty => ty, :x => x, :y => y}
end
# A DDA-direct search circle interpolator unrolled for each octant. Optimal and impure
def arc_unrolled(theta, angular_travel, radius)
radius = radius
x = (sin(theta)*radius).round
y = (cos(theta)*radius).round
angular_direction = angular_travel.sign
tx = (sin(theta+angular_travel)*(radius-0.5)).floor
ty = (cos(theta+angular_travel)*(radius-0.5)).floor
f = (x**2 + y**2 - radius**2).round
x2 = 2*x
y2 = 2*y
dx = (y==0) ? -x.sign : y.sign*angular_direction
dy = (x==0) ? -y.sign : -x.sign*angular_direction
max_steps = (angular_travel.abs*radius*2).floor
# Quandrants of the circls
# \ 1|2 /
# 8\ | / 3
# \|/
# ---------|-----------
# 7 /|\ 4
# / | \
# / 6 | 5 \
#
#
#
if angular_direction>0 # clockwise
if x.abs<y.abs # quad 1,2,6,5
if y>0 # quad 1,2
while x<0 # quad 1 x+,y+
x += 1
f += 1+x2
x2 += 2
f_diagonal = f + 1 + y2
if (f.abs >= f_diagonal.abs)
y += 1
y2 += 2
f = f_diagonal
end
end
while x>=0 # quad 2, x+, y-
x += 1
f += 1+x2
x2 += 2
f_diagonal = f + 1 - y2
if (f.abs >= f_diagonal.abs)
y -= 1
y2 -= 2
f = f_diagonal
end
end
end
if y<=0 # quad 6, 5
while x<0 # quad 6 x-, y+
x -= 1
f += 1-x2
x2 -= 2
f_diagonal = f + 1 + y2
if (f.abs >= f_diagonal.abs)
y += 1
y2 += 2
f = f_diagonal
end
end
while x>=0 # quad 5 x-, y-
x -= 1
f += 1-x2
x2 -= 2
f_diagonal = f + 1 - y2
if (f.abs >= f_diagonal.abs)
y -= 1
y2 -= 2
f = f_diagonal
end
end
end
# Quandrants of the circls
# \ 1|2 /
# 8\ | / 3
# \|/
# ---------|-----------
# 7 /|\ 4
# / | \
# / 6 | 5 \
else 3 # quad 3,4,7,8
if x>0 # quad 3,4
while y>0 # quad 3 x+,y+
x += 1
f += 1+x2
x2 += 2
f_diagonal = f + 1 + y2
if (f.abs >= f_diagonal.abs)
y += 1
y2 += 2
f = f_diagonal
end
end
while x>=0 # quad 2, x+, y-
x += 1
f += 1+x2
x2 += 2
f_diagonal = f + 1 - y2
if (f.abs >= f_diagonal.abs)
y -= 1
y2 -= 2
f = f_diagonal
end
end
end
if y<=0 # quad 6, 5
while x<0 # quad 6 x-, y+
x -= 1
f += 1-x2
x2 -= 2
f_diagonal = f + 1 + y2
if (f.abs >= f_diagonal.abs)
y += 1
y2 += 2
f = f_diagonal
end
end
while x>=0 # quad 5 x-, y-
x -= 1
f += 1-x2
x2 -= 2
f_diagonal = f + 1 - y2
if (f.abs >= f_diagonal.abs)
y -= 1
y2 -= 2
f = f_diagonal
end
end
end
else
y += dy
f += 1+y2*dy
y2 += 2*dy
f_diagonal = f + 1 + x2*dx
if (f.abs >= f_diagonal.abs)
x += dx
dy = -x.sign*angular_direction unless x == 0
x2 += 2*dx
f = f_diagonal
end
dx = y.sign*angular_direction unless y == 0
end
break if x*ty.sign*angular_direction>=tx*ty.sign*angular_direction && y*tx.sign*angular_direction<=ty*tx.sign*angular_direction
end
plot_pixel(tx+20, -ty+20, "o")
return {:tx => tx, :ty => ty, :x => x, :y => y}
end
end
test = CircleTest.new
test.init
#test.arc_clean(0, Math::PI*2, 5)
(1..10000).each do |r|
test.init
data = test.arc_supaoptimal(2.9, Math::PI*1, r)
if (data[:tx]-data[:x]).abs > 1 || (data[:ty]-data[:y]).abs > 1
puts "r=#{r} fails target control"
pp data
puts
end
end
# test.init
# data = test.arc_supaoptimal(1.1, -Math::PI, 19)
# pp data
#test.pure_arc(0,1,1,4)
test.show