Upload 6270 files

4083225 verified 6 months ago

3.85 kB

import {Adder} from "d3-array";
import {asin, atan2, cos, degrees, epsilon, epsilon2, hypot, radians, sin, sqrt} from "./math.js";
import noop from "./noop.js";
import stream from "./stream.js";

var W0, W1,
    X0, Y0, Z0,
    X1, Y1, Z1,
    X2, Y2, Z2,
    lambda00, phi00, // first point
    x0, y0, z0; // previous point

var centroidStream = {
  sphere: noop,
  point: centroidPoint,
  lineStart: centroidLineStart,
  lineEnd: centroidLineEnd,
  polygonStart: function() {
    centroidStream.lineStart = centroidRingStart;
    centroidStream.lineEnd = centroidRingEnd;
  },
  polygonEnd: function() {
    centroidStream.lineStart = centroidLineStart;
    centroidStream.lineEnd = centroidLineEnd;
  }
};

// Arithmetic mean of Cartesian vectors.
function centroidPoint(lambda, phi) {
  lambda *= radians, phi *= radians;
  var cosPhi = cos(phi);
  centroidPointCartesian(cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi));
}

function centroidPointCartesian(x, y, z) {
  ++W0;
  X0 += (x - X0) / W0;
  Y0 += (y - Y0) / W0;
  Z0 += (z - Z0) / W0;
}

function centroidLineStart() {
  centroidStream.point = centroidLinePointFirst;
}

function centroidLinePointFirst(lambda, phi) {
  lambda *= radians, phi *= radians;
  var cosPhi = cos(phi);
  x0 = cosPhi * cos(lambda);
  y0 = cosPhi * sin(lambda);
  z0 = sin(phi);
  centroidStream.point = centroidLinePoint;
  centroidPointCartesian(x0, y0, z0);
}

function centroidLinePoint(lambda, phi) {
  lambda *= radians, phi *= radians;
  var cosPhi = cos(phi),
      x = cosPhi * cos(lambda),
      y = cosPhi * sin(lambda),
      z = sin(phi),
      w = atan2(sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w), x0 * x + y0 * y + z0 * z);
  W1 += w;
  X1 += w * (x0 + (x0 = x));
  Y1 += w * (y0 + (y0 = y));
  Z1 += w * (z0 + (z0 = z));
  centroidPointCartesian(x0, y0, z0);
}

function centroidLineEnd() {
  centroidStream.point = centroidPoint;
}

// See J. E. Brock, The Inertia Tensor for a Spherical Triangle,
// J. Applied Mechanics 42, 239 (1975).
function centroidRingStart() {
  centroidStream.point = centroidRingPointFirst;
}

function centroidRingEnd() {
  centroidRingPoint(lambda00, phi00);
  centroidStream.point = centroidPoint;
}

function centroidRingPointFirst(lambda, phi) {
  lambda00 = lambda, phi00 = phi;
  lambda *= radians, phi *= radians;
  centroidStream.point = centroidRingPoint;
  var cosPhi = cos(phi);
  x0 = cosPhi * cos(lambda);
  y0 = cosPhi * sin(lambda);
  z0 = sin(phi);
  centroidPointCartesian(x0, y0, z0);
}

function centroidRingPoint(lambda, phi) {
  lambda *= radians, phi *= radians;
  var cosPhi = cos(phi),
      x = cosPhi * cos(lambda),
      y = cosPhi * sin(lambda),
      z = sin(phi),
      cx = y0 * z - z0 * y,
      cy = z0 * x - x0 * z,
      cz = x0 * y - y0 * x,
      m = hypot(cx, cy, cz),
      w = asin(m), // line weight = angle
      v = m && -w / m; // area weight multiplier
  X2.add(v * cx);
  Y2.add(v * cy);
  Z2.add(v * cz);
  W1 += w;
  X1 += w * (x0 + (x0 = x));
  Y1 += w * (y0 + (y0 = y));
  Z1 += w * (z0 + (z0 = z));
  centroidPointCartesian(x0, y0, z0);
}

export default function(object) {
  W0 = W1 =
  X0 = Y0 = Z0 =
  X1 = Y1 = Z1 = 0;
  X2 = new Adder();
  Y2 = new Adder();
  Z2 = new Adder();
  stream(object, centroidStream);

  var x = +X2,
      y = +Y2,
      z = +Z2,
      m = hypot(x, y, z);

  // If the area-weighted ccentroid is undefined, fall back to length-weighted ccentroid.
  if (m < epsilon2) {
    x = X1, y = Y1, z = Z1;
    // If the feature has zero length, fall back to arithmetic mean of point vectors.
    if (W1 < epsilon) x = X0, y = Y0, z = Z0;
    m = hypot(x, y, z);
    // If the feature still has an undefined ccentroid, then return.
    if (m < epsilon2) return [NaN, NaN];
  }

  return [atan2(y, x) * degrees, asin(z / m) * degrees];
}

	import {Adder} from "d3-array";
	import {asin, atan2, cos, degrees, epsilon, epsilon2, hypot, radians, sin, sqrt} from "./math.js";
	import noop from "./noop.js";
	import stream from "./stream.js";

	var W0, W1,
	X0, Y0, Z0,
	X1, Y1, Z1,
	X2, Y2, Z2,
	lambda00, phi00, // first point
	x0, y0, z0; // previous point

	var centroidStream = {
	sphere: noop,
	point: centroidPoint,
	lineStart: centroidLineStart,
	lineEnd: centroidLineEnd,
	polygonStart: function() {
	centroidStream.lineStart = centroidRingStart;
	centroidStream.lineEnd = centroidRingEnd;
	},
	polygonEnd: function() {
	centroidStream.lineStart = centroidLineStart;
	centroidStream.lineEnd = centroidLineEnd;
	}
	};

	// Arithmetic mean of Cartesian vectors.
	function centroidPoint(lambda, phi) {
	lambda = radians, phi = radians;
	var cosPhi = cos(phi);
	centroidPointCartesian(cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi));
	}

	function centroidPointCartesian(x, y, z) {
	++W0;
	X0 += (x - X0) / W0;
	Y0 += (y - Y0) / W0;
	Z0 += (z - Z0) / W0;
	}

	function centroidLineStart() {
	centroidStream.point = centroidLinePointFirst;
	}

	function centroidLinePointFirst(lambda, phi) {
	lambda = radians, phi = radians;
	var cosPhi = cos(phi);
	x0 = cosPhi * cos(lambda);
	y0 = cosPhi * sin(lambda);
	z0 = sin(phi);
	centroidStream.point = centroidLinePoint;
	centroidPointCartesian(x0, y0, z0);
	}

	function centroidLinePoint(lambda, phi) {
	lambda = radians, phi = radians;
	var cosPhi = cos(phi),
	x = cosPhi * cos(lambda),
	y = cosPhi * sin(lambda),
	z = sin(phi),
	w = atan2(sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w), x0 * x + y0 * y + z0 * z);
	W1 += w;
	X1 += w * (x0 + (x0 = x));
	Y1 += w * (y0 + (y0 = y));
	Z1 += w * (z0 + (z0 = z));
	centroidPointCartesian(x0, y0, z0);
	}

	function centroidLineEnd() {
	centroidStream.point = centroidPoint;
	}

	// See J. E. Brock, The Inertia Tensor for a Spherical Triangle,
	// J. Applied Mechanics 42, 239 (1975).
	function centroidRingStart() {
	centroidStream.point = centroidRingPointFirst;
	}

	function centroidRingEnd() {
	centroidRingPoint(lambda00, phi00);
	centroidStream.point = centroidPoint;
	}

	function centroidRingPointFirst(lambda, phi) {
	lambda00 = lambda, phi00 = phi;
	lambda = radians, phi = radians;
	centroidStream.point = centroidRingPoint;
	var cosPhi = cos(phi);
	x0 = cosPhi * cos(lambda);
	y0 = cosPhi * sin(lambda);
	z0 = sin(phi);
	centroidPointCartesian(x0, y0, z0);
	}

	function centroidRingPoint(lambda, phi) {
	lambda = radians, phi = radians;
	var cosPhi = cos(phi),
	x = cosPhi * cos(lambda),
	y = cosPhi * sin(lambda),
	z = sin(phi),
	cx = y0 * z - z0 * y,
	cy = z0 * x - x0 * z,
	cz = x0 * y - y0 * x,
	m = hypot(cx, cy, cz),
	w = asin(m), // line weight = angle
	v = m && -w / m; // area weight multiplier
	X2.add(v * cx);
	Y2.add(v * cy);
	Z2.add(v * cz);
	W1 += w;
	X1 += w * (x0 + (x0 = x));
	Y1 += w * (y0 + (y0 = y));
	Z1 += w * (z0 + (z0 = z));
	centroidPointCartesian(x0, y0, z0);
	}

	export default function(object) {
	W0 = W1 =
	X0 = Y0 = Z0 =
	X1 = Y1 = Z1 = 0;
	X2 = new Adder();
	Y2 = new Adder();
	Z2 = new Adder();
	stream(object, centroidStream);

	var x = +X2,
	y = +Y2,
	z = +Z2,
	m = hypot(x, y, z);

	// If the area-weighted ccentroid is undefined, fall back to length-weighted ccentroid.
	if (m < epsilon2) {
	x = X1, y = Y1, z = Z1;
	// If the feature has zero length, fall back to arithmetic mean of point vectors.
	if (W1 < epsilon) x = X0, y = Y0, z = Z0;
	m = hypot(x, y, z);
	// If the feature still has an undefined ccentroid, then return.
	if (m < epsilon2) return [NaN, NaN];
	}

	return [atan2(y, x) * degrees, asin(z / m) * degrees];
	}