Skip to content

Instantly share code, notes, and snippets.

@mathdoodle

mathdoodle/README.md

Last active September 23, 2016 19:54
Show Gist options
  • Save mathdoodle/eb1f40bd06787b53b8e0b364d498c9f1 to your computer and use it in GitHub Desktop.
Save mathdoodle/eb1f40bd06787b53b8e0b364d498c9f1 to your computer and use it in GitHub Desktop.
Download ZIP
orthoFramesToVersor algorithm
Raw
README.md

orthoFramesToVersor algorithm

Overview

Experiment in implementing this algorithm in TypeScript over the GeoCAS.Multivector type.

Raw
cube-fs.glsl
varying highp vec4 vColor;
void main(void) {
gl_FragColor = vColor;
}
Raw
cube-vs.glsl
attribute vec3 aPosition;
uniform vec3 uColor;
uniform float uOpacity;
uniform mat4 uModel;
uniform mat4 uProjection;
uniform mat4 uView;
varying highp vec4 vColor;
void main(void) {
gl_Position = uProjection * uView * uModel * vec4(aPosition, 1.0);
vColor = vec4(uColor, uOpacity);
}
Raw
displayError.ts
/**
* Displays an exception by writing it to a <pre> element.
*/
export default function displayError(e: any) {
const stderr = <HTMLPreElement>document.getElementById('error')
stderr.style.color = "#FF0000"
stderr.innerHTML = `${e}`
}
Raw
h3ga.ts
/**
* The Field Adapter for a number field.
*/
const fa = new GeoCAS.NumberFieldAdapter();
/**
* A (diagonal) metric for a 4D space.
* It doesn't matter whether we use +1 or -1 for the signature corresponding to e0
*/
const metric = [1, 1, 1, 1];
export const one = GeoCAS.getScalar(1, metric, fa);
export const e0 = GeoCAS.getBasisVector(0, metric, fa);
export const e1 = GeoCAS.getBasisVector(1, metric, fa);
export const e2 = GeoCAS.getBasisVector(2, metric, fa);
export const e3 = GeoCAS.getBasisVector(3, metric, fa);
const e0inverse = e0.inv();
/**
* Constructs the representation of a vector in Homogeneous 3D space.
* pos is the vector in 3D space.
* Returns the vector in 4D homogeneous space.
*/
export function vector(pos: EIGHT.VectorE3): GeoCAS.Multivector<number> {
return pos.x * e1 + pos.y * e2 + pos.z * e3;
}
export function toVectorE3(x: GeoCAS.Multivector<number>): EIGHT.VectorE3 {
return {x: x.scp(e1), y: x.scp(e2), z: x.scp(e3)};
}
/**
* Constructs the representation of a point in Homogeneous 3D space.
* pos is the vector in 3D space.
* Returns the point in 4D homogeneous space.
*/
export function point(pos: EIGHT.VectorE3): GeoCAS.Multivector<number> {
return e0 + vector(pos);
}
/**
* Translate an element of the Homogeneous model.
* X is the element of the model.
* t is the translation vector.
*/
export function translate(X: GeoCAS.Multivector<number>, t: GeoCAS.Multivector<number>): GeoCAS.Multivector<number> {
return X + t ^ (e0inverse << X);
}
/**
* Computes the string representation of a Homogeneous multivector.
*/
export function repr(m: GeoCAS.Multivector<number>): string {
return m.asString(['e0', 'e1', 'e2', 'e3']);
}
/**
* Computes the direction of a line element.
*/
export function direction(element: GeoCAS.Multivector<number>): EIGHT.VectorE3 {
// TODO: This probably could be generalized to compute the direction
// of points, lines, and planes by returning an EIGHT.Geometric3.
// The return types would be scalar, vector, bivector (respectively).
const a = e0inverse << element;
return toVectorE3(a);
}
/**
* The support is the point in the element (point, line, ...) that is closest to the origin.
*/
export function support(element: GeoCAS.Multivector<number>): EIGHT.VectorE3 {
// Versor inverse is supported so this function should cover most cases.
const d = (e0inverse << (e0 ^ element)).div(e0inverse << element);
return toVectorE3(d)
}
Raw
index.html
<!doctype html>
<html>
<head>
<!-- STYLES-MARKER -->
<style>
/* STYLE-MARKER */
</style>
<script src='https://jspm.io/[email protected]'></script>
<!-- SHADERS-MARKER -->
<!-- SCRIPTS-MARKER -->
</head>
<body>
<div id='container'>
<canvas id='canvas3D' width='500' height='500'></canvas>
<canvas id='canvas2D' width='500' height='500'></canvas>
</div>
<pre id='error'></pre>
<script>
// CODE-MARKER
</script>
<script>
System.import('./index.js')
</script>
</body>
</html>
Raw
index.ts
const zero = EIGHT.Geometric3.zero()
const e1 = EIGHT.Geometric3.e1()
const e2 = EIGHT.Geometric3.e2()
const e3 = EIGHT.Geometric3.e3()
import {point, repr, translate, vector} from './h3ga';
import {direction, support} from './h3ga';
// Comment out the following line to use the standard window.requestAnimationFrame
import requestAnimationFrame from './requestAnimationFrame';
import WireCube from './WireCube';
import WireLine from './WireLine';
const engine = new EIGHT.Engine('canvas3D')
.size(500, 500)
.clearColor(0.1, 0.1, 0.1, 1.0)
.enable(EIGHT.Capability.DEPTH_TEST)
const scene = new EIGHT.Scene(engine)
const ambients: EIGHT.Facet[] = []
const camera = new EIGHT.PerspectiveCamera();
//
// The following lines move the camera using the Homogeneous Model.
// It's a rather roundabout way of doing something quite simple, but the
// important point is that the translate function distinguishes points from vectors
// and also works for other elements in the model.
//
/**
* The position of the camera in homogeneous coordinates.
*/
const eye = translate(point(camera.eye), vector(1.5 * e3))
camera.eye.copyVector(support(eye))
ambients.push(camera)
const dirLight = new EIGHT.DirectionalLight()
ambients.push(dirLight)
const trackball = new EIGHT.TrackballControls(camera, window)
trackball.subscribe(engine.canvas)
trackball.noPan = true;
/**
* This graphic is just here to create a reference object.
*/
const wireCube = new WireCube(engine)
scene.add(wireCube)
/**
* sphere will represent points.
* We don't add it to the scene.
* Instead, we'll render manually in the animation loop.
*/
const sphere = new EIGHT.Sphere({engine})
sphere.radius = 0.03
/**
* wireLine will be used to visualize lines.
* We don't add it to the scene.
* Instead, we'll render manually in the animation loop.
*/
const wireLine = new WireLine(engine)
wireLine.color = EIGHT.Color.slateblue
const P = translate(point(+0.5 * e1), vector(zero))
const Q = translate(point(-0.5 * e1), vector(zero));
/**
* A line is the outer product of two points.
*/
const L = translate(P ^ Q, vector(zero));
const overlay = new EIGHT.Diagram3D('canvas2D', camera)
const animate = function(timestamp: number) {
engine.clear()
overlay.clear();
trackball.update()
dirLight.direction.copy(camera.look).sub(camera.eye)
wireLine.a.copyVector(direction(L))
wireLine.d.copyVector(support(L))
wireLine.render(ambients);
overlay.fillText('L = P ^ Q', wireLine.d)
sphere.X.copyVector(support(P))
sphere.color = EIGHT.Color.red
sphere.render(ambients)
overlay.fillText('P', sphere.X)
sphere.X.copyVector(support(Q))
sphere.color = EIGHT.Color.cobalt
sphere.render(ambients)
overlay.fillText('Q', sphere.X)
scene.draw(ambients)
requestAnimationFrame(animate)
}
requestAnimationFrame(animate)
Raw
line-fs.glsl
varying highp vec4 vColor;
void main(void) {
gl_FragColor = vColor;
}
Raw
line-vs.glsl
attribute float alpha;
uniform vec3 a;
uniform vec3 d;
uniform vec3 uColor;
uniform float uOpacity;
uniform mat4 uModel;
uniform mat4 uProjection;
uniform mat4 uView;
varying highp vec4 vColor;
void main(void) {
gl_Position = uProjection * uView * uModel * vec4(alpha * a + d, 1.0);
vColor = vec4(uColor, uOpacity);
}
Raw
math.ts
/**
* Universal exponential function for more natural notation.
*/
export function exp<T extends {exp: ()=>T}>(arg: T): T {
return arg.exp()
}
Raw
orthoFramesToVersor.ts
/**
* Create a type alias to make the algorithm a bit neater.
*/
type Multivector<T> = GeoCAS.Multivector<T>;
type Algebra<T> = GeoCAS.Algebra<T>;
const cos = GeoCAS.cosineOfAngleBetweenBlades;
/**
* Determines a versor for two frames related by an orthogonal transform.
*
* A and B are lists of corresponding vectors, representing the frames.
* A is the start frame, B is the final frame.
* V should be initialized to 1 in whatever multivector/field flavor you are using.
*
* This is an application of the Cartan-Dieudonne theorem that states:
* In n-D, any orthogonal transformation can be represented as at most n planar reflections.
*/
export default function orthoFramesToVersor<T>(A: Multivector<T>[], B: Multivector<T>[], vs: Multivector<T>[], algebra: Algebra<T>): Multivector<T>[] {
if (A.length > 0) {
const j = bestIndex(A, B, algebra);
const a = A[j];
const b = B[j];
const e = a.sub(b);
const field = algebra.field;
const eps = field.mulByNumber(field.one, 1e-6);
const cosMinusOne = cos(a, b).sub(algebra.one).scalarCoordinate();
const tooClose = field.lt(field.abs(cosMinusOne), eps);
if (tooClose) {
return orthoFramesToVersor(removeAt(A, j), removeAt(B, j), vs, algebra);
}
else {
const e2 = e.scp(e).scalarCoordinate();
const rvs = prepend(vs, e.divByScalar(algebra.field.sqrt(e2)));
// Don't let irrational number rounding errors from sqrt propagate...
return orthoFramesToVersor(removeAt(A, j).map(x => e.mul(x.mul(e)).neg().divByScalar(e2)), removeAt(B, j), rvs, algebra);
}
}
else {
return vs;
}
}
function prepend<T>(xs: T[], x: T): T[] {
const result: T[] = [];
result.push(x);
for (let i = 0; i < xs.length; i++) {
result.push(xs[i]);
}
return result;
}
/**
* Returns a copy of the array with the element at index removed.
*/
function removeAt<T>(xs: T[], index: number): T[] {
const result: T[] = [];
for (let i = 0; i < xs.length; i++) {
if (i !== index) {
result.push(xs[i]);
}
}
return result;
}
/**
* Determine the best vector for numerical stability.
*/
function bestIndex<T>(A: Multivector<T>[], B: Multivector<T>[], algebra: Algebra<T>): number {
const N = A.length;
let max = algebra.zero;
let idx = 0;
for (let k = 0; k < N; k++) {
const x = A[k].sub(B[k]);
const squaredNorm = x.scp(x.rev()).scalarCoordinate();
if (algebra.field.gt(squaredNorm, max.scalarCoordinate())) {
idx = k;
}
}
return idx;
}
Raw
package.json
{
"description": "orthoFramesToVersor algorithm",
"dependencies": {
"davinci-eight": "2.311.0",
"dat-gui": "0.5.0",
"geocas": "1.6.0"
},
"operatorOverloading": true,
"name": "copy-of-copy-of-copy-of-eight-starter-template",
"version": "0.1.0",
"keywords": [
"EIGHT",
"GeoCAS",
"Multivector",
"Homogeneous",
"h3ga",
"Cartan",
"Dieudonne"
],
"author": "David Geo Holmes"
}
Raw
requestAnimationFrame.ts
import displayError from './displayError'
/**
* Catches exceptions thrown in the animation callback and displays them.
* This function will have a slight performance impact owing to the try...catch statement.
* This function may be bypassed for production use by using window.requestAnimationFrame directly.
*/
export default function requestAnimationFrame(callback: FrameRequestCallback): number {
const wrapper: FrameRequestCallback = function(time: number) {
try {
callback(time)
}
catch(e) {
displayError(e)
}
}
return window.requestAnimationFrame(wrapper)
}
Raw
style.css
body { margin: 0; }
#container {
position: relative;
}
#canvas3D {
position: absolute;
left: 0px;
top: 0px;
width: 500px;
height: 500px;
}
#canvas2D {
position: absolute;
left: 0px;
top: 0px;
z-index: 10;
width: 500px;
height: 500px;
/**
* Allow events to go to the other elements
*/
pointer-events: none;
}
#stats { position: absolute; top: 0; left: 0; }
Raw
WireCube.ts
/**
* A Geometry for rendering a cube made from lines.
*/
class WireCubeGeometry implements EIGHT.Geometry {
private buffer: WebGLBuffer;
public data: Float32Array;
/**
*
*/
public invalid = true;
constructor(private engine: EIGHT.Engine) {
const gl = engine.gl;
const size = 1;
const L = size / 2;
this.data = new Float32Array([
-L, -L, -L, +L, -L, -L,
-L, +L, -L, +L, +L, -L,
-L, -L, +L, +L, -L, +L,
-L, +L, +L, +L, +L, +L,
-L, +L, +L, -L, +L, -L,
+L, +L, +L, +L, +L, -L,
-L, -L, +L, -L, -L, -L,
+L, -L, +L, +L, -L, -L,
-L, -L, -L, -L, +L, -L,
+L, -L, -L, +L, +L, -L,
-L, -L, +L, -L, +L, +L,
+L, -L, +L, +L, +L, +L
]);
this.buffer = gl.createBuffer();
}
bind(material: EIGHT.Material): void {
const gl = this.engine.gl;
gl.bindBuffer(gl.ARRAY_BUFFER, this.buffer);
}
unbind(material: EIGHT.Material): void {
const gl = this.engine.gl;
gl.bindBuffer(gl.ARRAY_BUFFER, null);
}
draw(material: EIGHT.Material): void {
const gl = this.engine.gl;
const aPosition = material.getAttribLocation('aPosition');
if (this.invalid) {
gl.bufferData(gl.ARRAY_BUFFER, this.data, EIGHT.Usage.STATIC_DRAW);
this.invalid = false;
}
gl.vertexAttribPointer(aPosition, 3, EIGHT.DataType.FLOAT, true, 0, 0);
gl.enableVertexAttribArray(aPosition);
gl.drawArrays(EIGHT.BeginMode.LINES, 0, 24);
gl.disableVertexAttribArray(aPosition);
}
}
export default class WireCube extends EIGHT.Mesh<WireCubeGeometry, EIGHT.Material> {
constructor(engine: EIGHT.Engine) {
super(new WireCubeGeometry(engine), new EIGHT.HTMLScriptsMaterial(['cube-vs', 'cube-fs'], document, [], engine), engine);
}
setUniforms() {
// We don't have any special uniforms so there is nothing to set.
// However, this is how you would do it...
const material = this.material;
// material.uniform3f('a', this.a.x, this.a.y, this.a.z);
material.release();
super.setUniforms();
// FIXME: This is just noise. It does not add any value.
return this;
}
}
Raw
WireLine.ts
/**
* A Geometry for rendering a parameterized line.
*/
class WireLineGeometry implements EIGHT.Geometry {
private buffer: WebGLBuffer;
public data: Float32Array;
/**
*
*/
public invalid = true;
constructor(private engine: EIGHT.Engine) {
const gl = engine.gl;
const size = 1;
this.data = new Float32Array([-1, 0, 0, +1]);
this.buffer = gl.createBuffer();
}
bind(material: EIGHT.Material): void {
const gl = this.engine.gl;
gl.bindBuffer(gl.ARRAY_BUFFER, this.buffer);
}
unbind(material: EIGHT.Material): void {
const gl = this.engine.gl;
gl.bindBuffer(gl.ARRAY_BUFFER, null);
}
draw(material: EIGHT.Material): void {
const gl = this.engine.gl;
const alpha = material.getAttribLocation('alpha');
if (this.invalid) {
gl.bufferData(gl.ARRAY_BUFFER, this.data, EIGHT.Usage.STATIC_DRAW);
this.invalid = false;
}
gl.vertexAttribPointer(alpha, 1, EIGHT.DataType.FLOAT, true, 0, 0);
gl.enableVertexAttribArray(alpha);
gl.drawArrays(EIGHT.BeginMode.LINES, 0, 4);
gl.disableVertexAttribArray(alpha);
}
}
export default class WireLine extends EIGHT.Mesh<WireLineGeometry, EIGHT.Material> {
/**
* The direction of the line.
*/
public a = EIGHT.Geometric3.e1().clone();
/**
* The support of the line.
*/
public d = EIGHT.Geometric3.zero().clone();
constructor(engine: EIGHT.Engine) {
super(new WireLineGeometry(engine), new EIGHT.HTMLScriptsMaterial(['line-vs', 'line-fs'], document, [], engine), engine);
}
setUniforms(): WireLine {
const material = this.material;
material.uniform3f('a', this.a.x, this.a.y, this.a.z);
material.uniform3f('d', this.d.x, this.d.y, this.d.z);
material.release();
super.setUniforms();
return this;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment