Unproject mouse to sphere met dome/pinhole pespectief

0.Unproject mouse to sphere.pde

Quaternion2 qNow2 = new Quaternion2();
Quaternion2 qTo2 = new Quaternion2();

float R = 600;
PShape shape;
PVector anchor3D = new PVector();

void settings() {
  fullScreen(P3D);
}

void setup() {
  sphereDetail(64);
  shape = createShape(SPHERE, R);
  shape.rotateY(HALF_PI);
  shape.setStroke(false);
  shape.setTexture(loadImage("earth.jpg"));
}

void update() {
  if (frameCount==2) {
    surface.setLocation(width/2-height/2, 0);
    surface.setSize(height, height);
  }
  qNow2 = Quaternion2.slerp(.1, qNow2, qTo2);
}

void draw() {
  update();

  background(10);
  lights();

  perspective(atan(0.5)*2, 1, 1200, 10000);
  camera(0, 0, -1200, 0, 0, 0, 0, 1, 0);
  scale(-1, 1, 1); // dome-flip

  stroke(128);
  noFill();
  PVector v = qNow2.getAxis();
  pushMatrix();
  rotate(qNow2.getAngle(), v.x, v.y, v.z);
  shape(shape);
  popMatrix();

  PVector hit = getMouseOnSphere(mouseX, mouseY);

  if (hit != null) {
    stroke(0, 255, 0);
    line(hit.x, hit.y, hit.z, anchor3D.x, anchor3D.y, anchor3D.z);
    pushMatrix();
    translate(hit.x, hit.y, hit.z);
    noStroke();
    fill(255, 0, 0);
    //sphereDetail(16);
    sphere(20);
    popMatrix();
  }
}

void mouseDragged() {
  PVector from = getMouseOnSphere(pmouseX, pmouseY);
  PVector to = getMouseOnSphere(mouseX, mouseY);
  drag(from, to, anchor3D, 1);
}

void mousePressed() {
  if (mouseButton==RIGHT) anchor3D = getMouseOnSphere(mouseX, mouseY);
}

void keyPressed() {
  if (key=='x') anchor3D.mult(0);
}

void drag(PVector _from, PVector _to, PVector _anchor, float scalar) { //-1..1
  if (_anchor==null || _from==null || _to==null) return;
  PVector from = _from.sub(_anchor).normalize();
  PVector to = _to.sub(_anchor).normalize();
  PVector axis = from.cross(to);
  float half_angle = asin(from.dot(to))/2;
  if (Float.isNaN(half_angle)) return;
  Quaternion2 q = new Quaternion2(cos(half_angle), axis.x * sin(half_angle), axis.y * sin(half_angle), axis.z * sin(half_angle));
  q = Quaternion2.slerp(scalar, new Quaternion2(1, 0, 0, 0), q);
  qTo2.mult(q);
  qTo2.normalize();
}

Quaternion2.pde

// static
public static class Quaternion2 implements PConstants {
  float W, X, Y, Z;

  Quaternion2() {
    set(1, 0, 0, 0);
  }

  Quaternion2(float w, float x, float y, float z) {
    set(w, x, y, z);
  }

  Quaternion2(float w, PVector axis) {
    set(w, axis.x, axis.y, axis.z);
  }

  Quaternion2 set(float w, PVector axis) {
    set(w, axis.x, axis.y, axis.z);
    return this;
  }

  void set(float w, float x, float y, float z) {
    W = w;
    X = x;
    Y = y;
    Z = z;
  }

  Quaternion2 mult(Quaternion2 q) {
    float x = q.W * X + q.X * W + q.Y * Z - q.Z * Y;
    float y = q.W * Y - q.X * Z + q.Y * W + q.Z * X;
    float z = q.W * Z + q.X * Y - q.Y * X + q.Z * W;
    float w = q.W * W - q.X * X - q.Y * Y - q.Z * Z;
    set(w, x, y, z);
    return this;
  }

  void multiplyScalar(float s) {
    set(W*s, X*s, Y*s, Z*s);
  }

  Quaternion2 copy() {
    return new Quaternion2(W, X, Y, Z);
  }

  PVector applyTo(PVector v) {
    // nVidia SDK implementation
    PVector uv, uuv;
    PVector qvec = new PVector(X, Y, Z); //_v.x, _v.y, _v.z);
    uv = qvec.cross(v);    //uv = qvec ^ v;
    uuv = qvec.cross(uv);    //uuv = qvec ^ uv;
    uv.mult(2.0f * W);
    uuv.mult(2.0f);
    v.add(uv);
    v.add(uuv);
    return v;
  }

  Quaternion2 normalize() {
    float norme = PApplet.sqrt(W*W + X*X + Y*Y + Z*Z);
    if (norme == 0.0f) {
      W = 1.0f;
      X = Y = Z = 0.0f;
    } else {
      float recip = 1.0f/norme;
      W *= recip;
      X *= recip;
      Y *= recip;
      Z *= recip;
    }
    return this;
  }

  float getAngle() {
    float sinhalfangle = PApplet.sqrt(X*X+Y*Y+Z*Z);
    return 2.0f * PApplet.atan2(sinhalfangle, W);
  }

  PVector getAxis() {
    float sinhalfangle = PApplet.sqrt(X*X+Y*Y+Z*Z);
    if (sinhalfangle>0) {
      PVector axis = new PVector(X, Y, Z);
      axis.div(sinhalfangle);
      return axis;
    } else return new PVector(0, 0, 1);
  }

  /// Set the elements of the Quat to represent a rotation of angle around the axis (x,y,z)
  static Quaternion2 fromRotate(float angle, float x, float y, float z ) {  // angle in Radians!
    float epsilon = 0.0000001f;
    float len = PApplet.sqrt( x * x + y * y + z * z );
    if (len < epsilon) return new Quaternion2(); // if ~zero length axis, so reset rotation to zero.
    float inversenorm  = 1.0f / len;
    float coshalfangle = PApplet.cos( 0.5f * angle );
    float sinhalfangle = PApplet.sin( 0.5f * angle );
    float _x = x * sinhalfangle * inversenorm;
    float _y = y * sinhalfangle * inversenorm;
    float _z = z * sinhalfangle * inversenorm;
    float _w = coshalfangle;
    return new Quaternion2(_w, _x, _y, _z); //FIXED
  }

  //adapted from https://github.com/mrdoob/three.js/blob/707b44a18c30161efdd2023247af88a4bde8d302/src/math/Quaternion.js#L350-L360
  static Quaternion2 fromUnitVectors(PVector from, PVector to) {
    PVector v1 = new PVector();
    float EPS = 0.000001f;
    float r = from.dot(to) + 1;
    if (r<EPS) {
      r = 0;
      if (PApplet.abs(from.x) > PApplet.abs(from.z)) v1.set(-from.y, from.x, 0);
      else v1.set(0, -from.z, from.y);
    } else {
      v1 = from.cross(to);
    }
    return new Quaternion2(r, v1).normalize();
  }

  static Quaternion2 fromVectors(PVector from, PVector to) {
    return fromUnitVectors(from.copy().normalize(), to.copy().normalize());
  }

  /// Spherical Linear Interpolation
  /// As t goes from 0 to 1, the Quat object goes from "from" to "to"
  /// Reference: Shoemake at SIGGRAPH 89
  /// See also
  /// http://www.gamasutra.com/features/programming/19980703/quaternions_01.htm
  static Quaternion2 slerp(float t, Quaternion2 from, Quaternion2 to) {
    float epsilon = 0.00001f;
    float omega, cosomega, sinomega, scale_from, scale_to ;

    Quaternion2 quatTo = to.copy();

    // this is a dot product
    cosomega = from.X*to.X + from.Y*to.Y + from.Z*to.Z + from.W*to.W;

    if ( cosomega < 0.0f ) {
      cosomega = -cosomega;
      quatTo.X *= -1;
      quatTo.Y *= -1;
      quatTo.Z *= -1;
      quatTo.W *= -1;
    }

    if ( (1.0f - cosomega) > epsilon ) {
      omega = PApplet.acos(cosomega) ; // 0 <= omega <= Pi (see man acos)
      sinomega = PApplet.sin(omega) ;  // this sinomega should always be +ve so
      // could try sinomega=sqrt(1-cosomega*cosomega) to avoid a sin()?
      scale_from = PApplet.sin((1.0f - t) * omega) / sinomega ;
      scale_to = PApplet.sin(t * omega) / sinomega ;
    } else {
      /* --------------------------------------------------
       The ends of the vectors are very close
       we can use simple linear interpolation - no need
       to worry about the "spherical" interpolation
       -------------------------------------------------- */
      scale_from = 1.0f - t ;
      scale_to = t ;
    }

    //add
    return new Quaternion2(
      from.W * scale_from + quatTo.W * scale_to,
      from.X * scale_from + quatTo.X * scale_to,
      from.Y * scale_from + quatTo.Y * scale_to,
      from.Z * scale_from + quatTo.Z * scale_to
      );
  }

  public String toString() {
    return W + "," + X + "," + Y + "," + Z;
  }

  PVector toEulerAngle() { //const Quaterniond& q, double& roll, double& pitch, double& yaw) {
    PVector rollPitchYaw = new PVector();

    //roll (x-axis rotation)
    float sinr_cosp = +2.0f * (W*X + Y*Z);
    float cosr_cosp = +1.0f - 2.0f * (X*X + Y*Y);
    rollPitchYaw.x = PApplet.atan2(sinr_cosp, cosr_cosp);

    // pitch (y-axis rotation)
    float sinp = +2.0f * (W*Y - Z*X);
    if (PApplet.abs(sinp) >= 1)
      rollPitchYaw.y = sinp<0 ? -HALF_PI : HALF_PI; //copysign(M_PI / 2, sinp); // use 90 degrees if out of range
    else
      rollPitchYaw.y = PApplet.asin(sinp);

    // yaw (z-axis rotation)
    float siny_cosp = +2.0f * (W*Z + X*Y);
    float cosy_cosp = +1.0f - 2.0f * (Y*Y + Z*Z);
    rollPitchYaw.z = PApplet.atan2(siny_cosp, cosy_cosp);

    return rollPitchYaw;
  }

  float getYaw() {
    float siny_cosp = +2.0f * (W*Z + X*Y);
    float cosy_cosp = +1.0f - 2.0f * (Y*Y + Z*Z);
    return PApplet.atan2(siny_cosp, cosy_cosp);
  }

  float getLength2() {
    return X*X + Y*Y + Z*Z + W*W;
  }

  Quaternion2 div(Quaternion2 q) {
    mult(q.getInverted()); ///new Quaternion().invert(q));
    return this;
  }

  Quaternion2 getInverted() {
    float dot = getLength2();
    if (dot!=0) dot = 1/dot;
    return new Quaternion2(W * dot, -X * dot, -Y * dot, -Z * dot);
  }

  PVector getAppliedTo(PVector p) {
    return applyTo(p.copy());
  }

  boolean isIdentity() {
    return W==1 && X==0 && Y==0 && Z==0;
  }

  float getHeading() {
    PVector v = new PVector(0, 0, 1);
    applyTo(v);
    return atan2(v.y, v.x);
  }

  void reset() {
    set(1, 0, 0, 0);
  }
}

Utils.pde

import processing.opengl.*;
import processing.opengl.PGraphics3D;

PVector getMouseOnSphere(float x, float y) { //0..1200
  PVector[] ray = getMouseRay(x, y);
  return intersectRaySphere(ray[0], ray[1], new PVector(0, 0, 0), R);
}

// bouw ray uit muis via inverse(projection * modelview)
PVector[] getMouseRay(float x, float y) {
  PGraphics3D pg = (PGraphics3D) g;
  PMatrix3D proj = pg.projection.get();
  PMatrix3D mv   = pg.modelview.get();
  PMatrix3D pmv = proj.get();
  pmv.apply(mv);
  pmv.invert();
  PVector near = unproject(x, y, -1, pmv); // near in NDC
  PVector far  = unproject(x, y, 1, pmv); // far in NDC
  return new PVector[]{ near, far };
}

// scherm → wereld, gegeven inverse(P * MV)
PVector unproject(float sx, float sy, float ndcZ, PMatrix3D invPMV) {
  float x = sx / (float)width  * 2.0f - 1.0f;
  float y = (height - sy) / (float)height * 2.0f - 1.0f;
  float z = ndcZ;
  float[] in  = { x, y, z, 1 };
  float[] out = new float[4];
  invPMV.mult(in, out);
  if (out[3] == 0) return null;
  return new PVector(out[0] / out[3], out[1] / out[3], out[2] / out[3]);
}

// standaard ray–sphere intersectie
PVector intersectRaySphere(PVector p0, PVector p1, PVector center, float radius) {
  PVector dir = PVector.sub(p1, p0);
  dir.normalize();

  PVector oc = PVector.sub(p0, center);

  float a = dir.dot(dir);
  float b = 2 * oc.dot(dir);
  float c = oc.dot(oc) - radius * radius;

  float disc = b*b - 4*a*c;
  if (disc < 0) return null;

  float sqrtDisc = sqrt(disc);
  float t1 = (-b - sqrtDisc) / (2*a);
  float t2 = (-b + sqrtDisc) / (2*a);

  float t = Float.MAX_VALUE;
  if (t1 > 0 && t1 < t) t = t1;
  if (t2 > 0 && t2 < t) t = t2;
  if (t == Float.MAX_VALUE) return null;

  return PVector.add(p0, PVector.mult(dir, t));
}