gaussdensity()

in packages/miew/src/gfx/geometries/QuickSurfGeometry.js [32:116]
69 lines of code
10 McCabe index (conditional complexity)

  gaussdensity(surface, packedArrays, atomicNum, params) {
    const numAtoms = packedArrays.posRad.length / 4;
    const { posRad, colors } = packedArrays;
    const { numVoxels } = this;
    const { radScale, gaussLim, gridSpacing } = params;
    const invIsoValue = 1.0 / params.isoValue;
    const invGridSpacing = 1.0 / gridSpacing;
    const maxVoxelX = numVoxels[0] - 1;
    const maxVoxelY = numVoxels[1] - 1;
    const maxVoxelZ = numVoxels[2] - 1;
    // TODO is densityMap and volTexMap initialized?

    const { volMap, volTexMap } = surface;
    const volData = volMap.getData();
    const strideX = volMap.getStrideX();

    const volTexData = volTexMap.getData();
    const texStrideX = volTexMap.getStrideX();

    let atomWeightData;
    if (this._visibilitySelector != null) {
      atomWeightData = surface.atomWeightMap.getData();
    }

    const { atomMap } = surface;

    for (let i = 0; i < numAtoms; ++i) {
      const ind = i * 4;
      const scaledRad = posRad[ind + 3] * radScale;
      const atomicNumFactor = atomicNum === null ? 1.0 : atomicNum[i];
      const radInv = 1 / (2 * scaledRad * scaledRad);
      let radLim = gaussLim * scaledRad;
      const radLim2 = radLim * radLim;
      radLim *= invGridSpacing;

      let tmp = posRad[ind] * invGridSpacing;
      const xMin = Math.max((tmp - radLim) | 0, 0);
      const xMax = Math.min((tmp + radLim) | 0, maxVoxelX);
      tmp = posRad[ind + 1] * invGridSpacing;
      const yMin = Math.max((tmp - radLim) | 0, 0);
      const yMax = Math.min((tmp + radLim) | 0, maxVoxelY);
      tmp = posRad[ind + 2] * invGridSpacing;
      const zMin = Math.max((tmp - radLim) | 0, 0);
      const zMax = Math.min((tmp + radLim) | 0, maxVoxelZ);

      let dz = zMin * gridSpacing - posRad[ind + 2];
      for (let z = zMin; z <= zMax; ++z, dz += gridSpacing) {
        let dy = yMin * gridSpacing - posRad[ind + 1];
        for (let y = yMin; y <= yMax; ++y, dy += gridSpacing) {
          const dy2dz2 = dy * dy + dz * dz;

          if (dy2dz2 >= radLim2) {
            continue;
          }

          let addr = volMap.getDirectIdx(xMin, y, z);
          let texAddr = volTexMap.getDirectIdx(xMin, y, z);
          let dx = xMin * gridSpacing - posRad[ind];
          for (let x = xMin; x <= xMax; ++x, dx += gridSpacing, addr += strideX, texAddr += texStrideX) {
            const r2 = dx * dx + dy2dz2;
            const expVal = -r2 * radInv;

            let density = Math.exp(expVal) * atomicNumFactor;

            // store most relevant atom (with highest density)
            if (this._visibilitySelector != null
              && density > atomWeightData[addr]) { // NOSONAR
              atomWeightData[addr] = density;
              // we use same index into atom map and atomWeightMap
              atomMap[addr] = packedArrays.atoms[i];
            }

            volData[addr] += density;

            // TODO check for volTexMap routine?
            density *= invIsoValue;
            const colInd = i * 3;
            volTexData[texAddr] += density * colors[colInd];
            volTexData[texAddr + 1] += density * colors[colInd + 1];
            volTexData[texAddr + 2] += density * colors[colInd + 2];
          }
        }
      }
    }
  }