Brownian Motion on 3D Sphere

Stixan

I'm trying to make a brownian motion with a fixed stepping size on the surface of a 3D sphere. I want to produce a result similar to the animations on this page.

Turns out to be a real headscratcher. Best bet so far is the equations found in the section "Destination point given distance and bearing from start point" on this rather excellent page. Implementing these, I notice that the random walker seems to get attracted towards the poles on the sphere, leaving the impression that this is not at all a random walk.

Any tips, tricks, links, explanations or feedback/improvements to my code (provided below) is highly appreciated!

Crawler c;

void setup() {
  size(800, 800, P3D);
  c = new Crawler(200.0, 10.0);
}

void draw() {
  background(0);
  translate(width/2, height/2, 0);
  rotateY(frameCount * 0.01);
  c.display();
  c.update();
}

class Crawler {

  ArrayList<Point3D> positions;
  float radius, distance;
  float theta, phi;
  float bearing;

  Crawler(float _radius, float _stepSize) {
    radius = _radius;
    distance = _stepSize/radius;
    positions = new ArrayList<Point3D>();
    theta = random(-PI, PI);
    phi = random(-PI, PI);
    bearing = random(-PI, PI);
  }

  void update() {
    // Calculate new theta and phi values
    float phi2 = asin(sin(phi)*cos(distance)+cos(phi)*sin(distance)*cos(bearing));
    float theta2 = theta+atan2(sin(bearing)*sin(distance)*cos(phi), cos(distance)-sin(phi)*sin(phi2));

    // Calculate new bearing on destination
    float deltaTheta = theta-theta2;
    float tempy = sin(deltaTheta)*cos(phi);
    float tempx = cos(phi2)*sin(phi)-sin(phi2)*cos(phi)*cos(deltaTheta);
    float delta = atan2(tempy, tempx);
    float finalBearing = delta+PI;

    // Checking that the theta value are within the -PI to PI range
    if (theta2 < -PI) theta2 += TWO_PI;
    if (theta2 > PI) theta2-= TWO_PI;

    // Convert lat/lon to x/y/z positions and push position into array
    float x = radius * sin(theta2) * cos(phi2);
    float y = radius * sin(theta2) * sin(phi2);
    float z = radius * cos(theta2);
    positions.add(new Point3D(x, y, z));

    // Update theta and phi with new values
    theta = theta2;
    phi = phi2;

    // Update bearing with new value, add a random course and keep within 0 to TWO_PI range
    bearing = finalBearing + random(-HALF_PI, HALF_PI); // random(TWO_PI) = any direction possible
    if (bearing<0) bearing += TWO_PI;
    if (bearing>TWO_PI) bearing -= TWO_PI;
  }

  void display() {
    stroke(255);
    noFill();
    int counter = 0;
    float prevx = 0, prevy = 0, prevz = 0;
    for (Point3D p : positions) {
      if (counter > 1) {
        line(prevx, prevy, prevz, p.x, p.y, p.z);
      }
      prevx = p.x;
      prevy = p.y;
      prevz = p.z;
      counter++;
    }
  }

  class Point3D {
    float x, y, z;
    Point3D(float _x, float _y, float _z) {
      x = _x;
      y = _y;
      z = _z;
    }
  }

}

quark

Looks really nice.

I notice that the random walker seems to get attracted towards the poles on the sphere, leaving the impression that this is not at all a random walk.

The angular distance between points in the east-west position at the poles will be smaller that that at the equator for the same change in longitude. So it will take longer for a walker to leave a pole. You could overcome this using a torus instead of a sphere.

Rather than a random walk you might try the wander steering behaviour

Stixan

Qualified input. Thank you quark.

While I understand the difference between a torus and a sphere, I'm not sure how this would solve the problem points clumping together at certain parts of the torus and still cause "polar drift". Have you got any examples that involve using a torus instead of a sphere?

I'll look into the wander behaviour as you suggest.

Chrisir

point 3D is PVector in processing (see reference) btw.

quark

Imagine a 2D grid of regularly spaced points and you wrap it round a sphere. The points near the poles are closer together than at the equater, and at the poles they are touching each other.

If you wrap the grid of the toruus the points on the inside of the ring are closer together than the ones on the outside but to a lesser degree provided the ring radius is significantly larger than the tube radius.

This is a simplified version of your code but modified to work on a torus.

Crawler c;

void setup() {
  size(800, 800, P3D);
  c = new Crawler(180, 60);
}

void draw() {
  background(0);
  translate(width/2, height/2, 0);
  rotateX(frameCount * 0.01);
  rotateZ(frameCount * 0.002);
  c.display();
  c.update();
}

class Crawler {

  ArrayList<Point3D> positions;
  float tubeRad, ringRad;
  float theta, phi;
  float deltaPhi, deltaTheta;

  Crawler(float _ringRad, float _tubeRad) {
    tubeRad = _tubeRad;
    ringRad = _ringRad;
    positions = new ArrayList<Point3D>();
    phi = random(-PI, PI);
    theta = random(-PI, PI);
    deltaPhi = PI/15; // 12 degree
    deltaTheta = PI/30; // 6 degree
  }

  float fixAngle(float angle) {
    while (angle < -PI) angle += TWO_PI;
    while (angle > PI) angle -= TWO_PI;
    return angle;
  }

  void update() {
    // Calculate new theta and phi values
    phi += random(-deltaPhi, deltaPhi);
    theta += random(-deltaTheta, deltaTheta);
    phi = fixAngle(phi);
    theta = fixAngle(theta);

    // Convert lat/lon to x/y/z positions and push position into array
    float temp = ringRad + tubeRad * cos(phi);
    float px = temp * cos(theta);
    float py = tubeRad * sin(phi);
    float pz = temp * sin(theta);
    positions.add(new Point3D(px, py, pz));
  }

  void display() {
    stroke(255);
    noFill();
    int counter = 0;
    float prevx = 0, prevy = 0, prevz = 0;
    for (Point3D p : positions) {
      if (counter > 1) {
        line(prevx, prevy, prevz, p.x, p.y, p.z);
      }
      prevx = p.x;
      prevy = p.y;
      prevz = p.z;
      counter++;
    }
  }

  class Point3D {
    float x, y, z;
    Point3D(float _x, float _y, float _z) {
      x = _x;
      y = _y;
      z = _z;
    }
  }
}

Chrisir

can't we compensate the pole effect by modyfing a factor for random?

When closer to pole, less factor?

Stixan

Appreciate you taking the time to modify my example, quark. Very useful. Only reason I've chosen to reject this as an answer, is because my initial question was about a sphere not a torus. If, by some clever topological math, you can transform the torus shape into a perfect sphere and still retain the even spread in the crawlers path, I'll be chuffed to bits.

Chrisir, you might be on to something here. I've thought about this myself, but as soon as I start to fiddle with the formulas, I very quickly screw things up. At least when using a combination of spherical/cartesian coordinates. That's why I'm considering switching to using vectors, but that's a whole new ballgame to me. Found this example by Daniel Shiffman that I'm currently trying to wrangle into doing what I've described above.

G_D

Not sure if this will help, but I think your question essentially boils down to placing equidistant points on a sphere (as points otherwise get closer at each pole as quark has already pointed out). Anyway, from what I can tell, this is quite hard to achieve (in a purely mathematical sense) but it seems there are ways to approximate this. See code below...at least to me this seems to give a more even distribution of points around a sphere. Perhaps you can adapt to your crawler?

float radius = 200.0;
float x, y, z;

void setup() {
size(500, 500, P3D);
sphereDetail(30);
}

void draw() {
  background(50);
  stroke(255);
  translate(width/2, height/2);
  rotateX(radians(mouseY));
  rotateY(radians(mouseX));
  noFill();
  stroke(255,0,0);
  sphere(195); // just to show clustering at poles...

   for(float phi = 0.0; phi <= PI; phi += PI/30) {
    for(float theta = 0.0; theta <= TWO_PI; theta += 1/(15*sin(phi))) {

 //adapted from https://www.physicsforums.com/threads/uniform-grid-of-points-on-a-sphere.649761/

      x = radius * sin(phi) * cos(theta);
      y = -radius * cos(phi);
      z = radius * sin(phi) * sin(theta);

      pushMatrix();
      translate(x, y, z);
      noStroke();
      fill(255);
      box(3);
      popMatrix();
    }

  }

}

Stanlepunk

I work around the projection sphere in Polar and Cartesian mode maybe this piece of code can help you. You can remplace the class Vec3 by PVector or downlod the full code here, if you want : https://github.com/StanLepunK/Math

 /**
SPHERE PROJECTION
*/
/**
// FIBONACCI SPHERE PROJECTION CARTESIAN
*/
Vec3 [] list_cartesian_fibonacci_sphere (int num, float step, float root) {
  float root_sphere = root *num ;
  Vec3 [] list_points = new Vec3[num] ;
  for (int i = 0; i < list_points.length ; i++) list_points[i] = distribution_cartesian_fibonacci_sphere(i, num, step, root_sphere) ;
  return list_points ;
}
/*
float root_cartesian_fibonnaci(int num) {
  return random(1) *num ;
}
*/

Vec3 distribution_cartesian_fibonacci_sphere(int n, int num, float step, float root) {
  if(n<num) {
    float offset = 2. / num ;
    float y  = (n *offset) -1 + (offset / 2.) ;
    float r = sqrt(1 - pow(y,2)) ;
    float phi = ((n +root)%num) * step ;

    float x = cos(phi) *r ;
    float z = sin(phi) *r ;

    return Vec3(x,y,z) ;
  } else return Vec3() ;
}

/**
// POLAR PROJECTION FIBONACCI SPHERE
*/
Vec2 [] list_polar_fibonacci_sphere(int num, float step) {
  Vec2 [] list_points = new Vec2[num] ;
  for (int i = 0; i < list_points.length ; i++) list_points[i] = distribution_polar_fibonacci_sphere(i, num, step) ;
  return list_points ;
}
Vec2 distribution_polar_fibonacci_sphere(int n, int num, float step) {
  if(n<num) {
    float longitude = PHI *TAU *n;
    longitude /= step ;
    // like a normalization of the result ?
    longitude -= floor(longitude); 
    longitude *= TAU;
    if (longitude > PI)  longitude -= TAU;
    // Convert dome height (which is proportional to surface area) to latitude
    float latitude = asin(-1 + 2*n/(float)num);
    return Vec2(longitude, latitude) ;
  } else return Vec2() ;

}

Stixan

Thank you G_D! The formula you found on physicsforum.com was exactly what I needed. For anyone that might be interested, I'm sharing a sloppy yet working version my code that solves my initial question. Thank you all for your input. Internet high fives!

ArrayList<Crawler> crawlers;
ArrayList<Point3D> points;

int detail;
float radius, distance;

void setup() {
  size(500, 500, P3D);
  randomize();
}

void draw() {
  background(60);
  translate(width/2, height/2);
  rotateY(frameCount * 0.01);
  // Display crawlers
  for (Crawler c : crawlers) {
    c.update();
    c.display();
  }
}

void keyPressed() {
  if (key == ' ') randomize();
}

void randomize() {
  // Randomize sphere values
  detail = (int)random(10, 60);
  radius = random(100, 200.0);
  distance = random(TWO_PI / detail * radius, radius);
  // Calculate sphere points
  calcPoints(radius, detail);
  // Add a random number of crawlers
  crawlers = new ArrayList<Crawler>();
  int n = (int)random(1, 20);
  for (int i = 0; i < n; i++) {
    crawlers.add(new Crawler());
  }
}

void calcPoints(float _radius, int _detail) {
  points = new ArrayList<Point3D>();
  // Precalculate position of points evenly distributed around sphere
  float x, y, z;
  for (float phi = 0.0; phi <= PI; phi += PI/_detail) {
    for (float theta = 0.0; theta <= TWO_PI; theta += 1/((_detail/2)*sin(phi))) {
      x = _radius * sin(phi) * cos(theta);
      y = -_radius * cos(phi);
      z = _radius * sin(phi) * sin(theta);
      points.add(new Point3D(x, y, z));
    }
  }
}

class Crawler {
  ArrayList<Point3D> path;
  int pathLength, currentPoint;
  Crawler() {
    path = new ArrayList<Point3D>();
    int currentPoint = (int)random(points.size());
    path.add(points.get(currentPoint));
    pathLength = (int)random(10, 200);
  }
  void update() {
    // Calculate array of points within allowed distance (NEEDS OPTIMIZING!)
    Point3D cp = points.get(currentPoint);
    ArrayList<Integer> jumpList = new ArrayList<Integer>();
    for (int j=0; j<points.size(); j++) {
      Point3D p = points.get(j);
      if (dist(cp.x, cp.y, cp.z, p.x, p.y, p.z) < distance) jumpList.add(j);
    }
    // Choose a random valid point and add its coordinates to the path
    if (jumpList.size() > 0) {
      int n = (int)random(jumpList.size());
      currentPoint = jumpList.get(n);
      Point3D np = points.get(currentPoint);
      if (path.size() == pathLength) path.remove(0);
      path.add(np);
    }
  }
  void display() {
    // Display path
    stroke(255, 200);
    noFill();
    beginShape();
    for (Point3D p : path) {
      curveVertex(p.x, p.y, p.z);
    }
    endShape();
  }
}

class Point3D {
  float x, y, z;
  Point3D(float _x, float _y, float _z) {
    x = _x;
    y = _y;
    z = _z;
  }
}

Chrisir

Well done!

rbrauer

You could do it with a PMatrix and its axis / angle rotate method too. Basically start with a point offset from the center by radius. Apply the rotation and for each additional point use the rotated position of the previous point as the next point.

PMatrix3D pm;
PVector pv;
ArrayList<PVector> pts;

void setup() {
  size(500,500,P3D);
  stroke(255);
  noFill();

  pm = new PMatrix3D();
  pv = new PVector(0,0,200);
  pts = new ArrayList<PVector>();
}

void draw() {
  background(0);
  translate(width*0.5, height*0.5);
  rotateY(frameCount * 0.01);

  PVector axis = PVector.random3D();
  pm.rotate(0.7, axis.x, axis.y, axis.z);
  PVector next = pm.mult(pv, null);
  pm.reset();
  pts.add(next);
  if(pts.size() > 50) pts.remove(0);
  pv = next.get();

  beginShape();
  for(PVector cur : pts)
    curveVertex(cur.x, cur.y, cur.z);
  endShape();
}

Stixan

Wow, rbrauer, that's truly impressive! It never occurred to me that I could use PMatrix3D. I'll be spending the next couple of hours dissecting your code to try and understand how it works. Tip of the hat to you. And thank you for sharing.