matrix math circumcircle

Author: monkstone

processing_app/library/vecmath/vec2d/circles.rb

@@ -0,0 +1,50 @@
+#!/usr/bin/env jruby
+require 'propane'
+
+class Circles < Propane::App
+  load_library :circles
+
+  def settings
+    size(800, 600, P2D)
+  end
+
+  ## To be overriden by the Presentation Code.
+  def setup
+    sketch_title 'Circles'
+    color_mode(HSB, 360, 100, 100, 100)
+    @c = rand(360)
+    @points = TPoints.new
+    3.times { @points << TPoint.new(Vec2D.new(rand(5..width - 5), rand(5..height - 5))) }
+    background 0
+    ellipse_mode(RADIUS)
+  end
+
+  def draw
+    fill(0, 0, 0)
+    no_stroke
+    rect(0, 0, width, height) if (frame_count % 8_000).zero?
+    @points.each do |point|
+      # change direction sometimes
+      point.direction Vec2D.random if rand > 0.96
+      point.update
+    end
+
+    # set the style of the circle
+    @dc = map1d(millis, 0..150_000, 0..360) # slowly changes hue
+    stroke((@c + @dc) % 360, 50, 100, 5)
+    no_fill
+
+    ## verifies if there is a circle and draw it
+    draw_circle @points unless @points.collinear?
+  end
+
+  def draw_circle(pts)
+    circumcircle = Circumcircle.new(@points.vec)
+    circumcircle.calculate
+    center_point = circumcircle.center
+    radius = circumcircle.radius
+    ellipse(center_point.x, center_point.y, radius, radius)
+  end
+end
+
+Circles.new

processing_app/library/vecmath/vec2d/circumcircle.md

@@ -0,0 +1,60 @@
+### Dealing with processing coordinate system ###
+
+For library see `library/circles/circles.rb`
+
+If you want to create Math Sketches in processing, you need to deal with peculiar coordinate systems, where the Y-axis is inverted in theory you should be able to do.
+
+```ruby
+scale(1, -1)
+translate(0, -height)
+```
+But that will mess up any text (plus you probably need to `push_matrix` and `pop_matrix`) so it is probably simpler to create a parallel coordinate system for the math, and translate that back to the screen (using the processing `map` function or in `propane` and `JRubyArt` use `map1d`).
+
+We have done this in `circumcircle_sketch.rb` or just accept the processing coordinate system as we have with `basic_cirmcumcircle_sketch.rb` (it is much simpler).
+
+### PVector limitations ###
+
+PVector is a 3D vector, that is often used as a 2D vector, to my mind that is just plain wrong. If you evaluate the cross product of a 2D vector you get a float (you may see somewhere that a cross product does not exist for 2D vectors, but it can be useful).  The cross product of PVector yields another PVector, so cannot be used in the calculation of the area of the triangle as defined by two vectors _cf_ Vec2D:-
+
+```ruby
+a = Vec2D.new(100, 0)
+b = Vec2D.new(0, 100)
+
+a.cross(b).abs == 10_000 # or twice the area of the triangle enclosed by a, b
+
+```
+Further we can use the cross product in a test for collinearity
+
+```ruby
+
+# given 3 points in 2D space
+a = Vec2D.new(0, 0)
+b = Vec2D.new(100, 100)
+c = Vec2D.new(200, 200)
+
+(a - b).cross(b - c) == 0 # the area of the triangle is zero, so a, b, c are collinear
+
+```
+
+Also because we were able to separate the logic, we were able to confidently re-factor the Barbara Almeida sketch to use Matrix math to determine the circumcenter
+
+### Matrix Math ###
+
+For detailed workings see [Circumcircle at Mathworld Wolfram.com][circumcircle]
+
+
+a = {{x<sub>1</sub> y<sub>1</sub> 1}, {x<sub>2</sub> y<sub>2</sub> 1}, {x<sub>3</sub> y<sub>3</sub> 1}}
+
+bx = -{{x<sub>1</sub><sup>2</sup> + y<sub>1</sub><sup>2</sup> y<sub>1</sub> 1}, {x<sub>2</sub><sup>2</sup> + y<sub>2</sub><sup>2</sup> y<sub>2</sub> 1}, {x<sub>3</sub><sup>2</sup> + y<sub>3</sub><sup>2</sup> y<sub>3</sub> 1}}
+
+by = {{x<sub>1</sub><sup>2</sup> + y<sub>1</sub><sup>2</sup> x<sub>1</sub> 1}, {x<sub>2</sub><sup>2</sup> + y<sub>2</sub><sup>2</sup> x<sub>2</sub> 1}, {x<sub>3</sub><sup>2</sup> + y<sub>3</sub><sup>2</sup> x<sub>3</sub> 1}}
+
+xo = -bx / 2 * a
+
+yo = -by / 2 * a
+
+
+[circumcircle]:http://mathworld.wolfram.com/Circumcircle.html
+
+
+

processing_app/library/vecmath/vec2d/library/circles/circles.rb

@@ -0,0 +1,3 @@
+require_relative 'lib/circumcircle'
+require_relative 'lib/t_points'
+require_relative 'lib/triangle_point'

processing_app/library/vecmath/vec2d/library/circles/lib/circumcircle.rb

@@ -0,0 +1,43 @@
+# Circumcircle from 3 points
+require 'matrix'
+
+class Circumcircle
+  attr_reader :center, :radius, :points
+  def initialize(points)
+    @points = points
+  end
+
+  def calculate
+    @center = Vec2D.new(-(bx / am), -(by / am))
+    @radius = center.dist(points[2]) # any point would do
+  end
+
+  private
+
+  # Matrix math see matrix_math.md and in detail
+  # http://mathworld.wolfram.com/Circumcircle.html
+
+  def am
+    2 * Matrix[
+      [points[0].x, points[0].y, 1],
+      [points[1].x, points[1].y, 1],
+      [points[2].x, points[2].y, 1]
+    ].determinant
+  end
+
+  def bx
+    -Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].y, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].y, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].y, 1]
+    ].determinant
+  end
+
+  def by
+    Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].x, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].x, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].x, 1]
+    ].determinant
+  end
+end

processing_app/library/vecmath/vec2d/library/circles/lib/t_points.rb

@@ -0,0 +1,33 @@
+# frozen_string_literal: true
+require 'forwardable'
+MAX_POINT = 3
+# A collection of a maximum of 3 points in the processing world
+# includes a collinearity test using Vec2D
+class TPoints
+  extend Forwardable
+  def_delegators(:@points, :each, :map, :size, :shift, :clear, :[])
+  include Enumerable
+
+  attr_reader :points
+
+  def initialize
+    @points = []
+  end
+
+  def <<(pt)
+    points << pt
+    shift if size > MAX_POINT
+  end
+
+  def collinear?
+    full? ? (vec[0] - vec[1]).cross(vec[1] - vec[2]).zero? : false
+  end
+
+  def vec
+    points.map { |point| point.pos }
+  end
+
+  def full?
+    points.length == MAX_POINT
+  end
+end

processing_app/library/vecmath/vec2d/library/circles/lib/triangle_point.rb

@@ -0,0 +1,43 @@
+# frozen_string_literal: true
+# particle and triangle point
+class TPoint
+  include Propane::Proxy
+  attr_reader :pos, :vel, :accel, :xbound, :ybound
+
+  def initialize(position)
+    @pos = position
+    @vel = Vec2D.new
+    @accel = Vec2D.random
+    @xbound = Boundary.new(0, width)
+    @ybound = Boundary.new(0, height)
+  end
+
+  def direction(acc)
+    # magnitude of the acceleration is proportional to the angle between
+    # acceleration and velocity
+    dif = accel.angle_between(vel)
+    dif = map1d(dif, 0..PI, 0.1..0.001)
+    @accel = acc * dif
+  end
+
+  def update
+    @vel += accel
+    @vel.set_mag(1.5) { vel.mag > 1.5 }
+    @pos += vel
+    check_bounds
+  end
+
+  private
+
+  def check_bounds
+    @vel.x *= -1 if xbound.exclude? pos.x
+    @vel.y *= -1 if ybound.exclude? pos.y
+  end
+end
+
+# we are looking for excluded values
+Boundary = Struct.new(:lower, :upper) do
+  def exclude?(val)
+    true unless (lower...upper).cover? val
+  end
+end