_forEachVoxelWithinRadius()

in packages/miew/src/chem/VoxelWorld.js [165:228]


  _forEachVoxelWithinRadius(center, radius, process) {
    const xRange = VoxelWorld._xRange;
    const yRange = VoxelWorld._yRange;
    const zRange = VoxelWorld._zRange;

    // switch to a faster method unless cell size is much smaller than sphere radius
    if (radius / this._cellInnerR < 10) {
      this._forEachVoxelWithinRadiusSimple(center, radius, process);
      return;
    }

    let rRangeXY;
    let rRangeX;
    let xVal;
    let yVal;
    let zVal;
    let isInsideX;
    let isInsideY;
    let isInsideZ;

    zRange.set(center.z - radius, center.z + radius);
    zRange.subScalar(this._box.min.z)
      .divideScalar(this._cellSize.z)
      .floor()
      .clampScalar(0, this._count.z - 1);

    for (let z = zRange.x; z <= zRange.y; ++z) {
      zVal = [this._box.min.z + z * this._cellSize.z,
        this._box.min.z + (z + 1) * this._cellSize.z];

      isInsideZ = (center.z - radius <= zVal[0]) && (zVal[1] <= center.z + radius);

      rRangeXY = _getSphereSliceRadiusRange(center, radius, zVal[0], zVal[1]);

      yRange.set(center.y - rRangeXY[1], center.y + rRangeXY[1]);
      yRange.subScalar(this._box.min.y)
        .divideScalar(this._cellSize.y)
        .floor()
        .clampScalar(0, this._count.y - 1);

      for (let y = yRange.x; y <= yRange.y; ++y) {
        yVal = [this._box.min.y + y * this._cellSize.y,
          this._box.min.y + (y + 1) * this._cellSize.y];

        isInsideY = (center.y - rRangeXY[0] <= yVal[0]) && (yVal[1] <= center.y + rRangeXY[0]);

        rRangeX = _getCircleSliceRadiusRange(center, rRangeXY[1], yVal[0], yVal[1]);

        xRange.set(center.x - rRangeX[1], center.x + rRangeX[1]);
        xRange.subScalar(this._box.min.x)
          .divideScalar(this._cellSize.x)
          .floor()
          .clampScalar(0, this._count.x - 1);

        for (let { x } = xRange; x <= xRange.y; ++x) {
          xVal = [this._box.min.x + x * this._cellSize.x,
            this._box.min.x + (x + 1) * this._cellSize.x];
          isInsideX = (center.x - rRangeX[0] <= xVal[0]) && (xVal[1] <= center.x + rRangeX[0]);

          process(x + this._count.x * (y + this._count.y * z), isInsideX && isInsideY && isInsideZ);
        }
      }
    }
  }