in packages/miew/src/gfx/geometries/IsoSurface.js [445:622]
setColorVolTex(colorMap, atomMap, atomWeightMap, visibilitySelector) {
let i;
let idx;
const numVerts = this._position.length / 3;
const vertices = this._position;
const origin = this._origin;
const dim = this._volumetricData.getDimensions();
const xs = dim[0] - 1;
const ys = dim[1] - 1;
const zs = dim[2] - 1;
const colorData = colorMap.getData();
const strideX = colorMap.getStrideX();
const strideY = colorMap.getStrideY();
const strideZ = colorMap.getStrideZ();
let atomWeightData;
let atomStrideX;
let atomStrideY;
let atomStrideZ;
if (visibilitySelector !== null) {
atomWeightData = atomWeightMap.getData();
atomStrideX = atomWeightMap.getStrideX();
atomStrideY = atomWeightMap.getStrideY();
atomStrideZ = atomWeightMap.getStrideZ();
}
const xInv = 1.0 / this._xAxis.x;
const yInv = 1.0 / this._yAxis.y;
const zInv = 1.0 / this._zAxis.z;
let atomLookup = [];
let atomWeights = [];
const colors = utils.allocateTyped(Float32Array, numVerts * 3);
function interp(mu, idx1, idx2, c) {
c[0] = (1 - mu) * colorData[idx1] + mu * colorData[idx2];
c[1] = (1 - mu) * colorData[idx1 + 1] + mu * colorData[idx2 + 1];
c[2] = (1 - mu) * colorData[idx1 + 2] + mu * colorData[idx2 + 2];
}
function collectWeight(ai, coefX, coefY, coefZ) {
const a = atomMap[ai]; // atomWeightMap is a scalar field, so index into atom map should be the same
if (a != null) {
atomLookup[a.index] = a;
const w = coefX * coefY * coefZ * atomWeightData[ai];
if (typeof atomWeights[a.index] === 'undefined') {
atomWeights[a.index] = w;
} else {
atomWeights[a.index] += w;
}
}
}
const vMap = utils.allocateTyped(Int32Array, numVerts);
let newVerCount = 0;
for (i = 0; i < numVerts; i++) {
const ind = i * 3;
const vx = (vertices[ind] - origin.x) * xInv;
const vy = (vertices[ind + 1] - origin.y) * yInv;
const vz = (vertices[ind + 2] - origin.z) * zInv;
const x = Math.min(Math.max(vx, 0), xs) | 0;
const y = Math.min(Math.max(vy, 0), ys) | 0;
const z = Math.min(Math.max(vz, 0), zs) | 0;
const mux = (vx - x);
const muy = (vy - y);
const muz = (vz - z);
if (visibilitySelector != null) {
// collect atom weights
atomLookup = [];
atomWeights = [];
idx = atomWeightMap.getDirectIdx(x, y, z);
collectWeight(idx, 1 - mux, 1 - muy, 1 - muz);
collectWeight(idx + atomStrideX, mux, 1 - muy, 1 - muz);
collectWeight(idx + atomStrideY, 1 - mux, muy, 1 - muz);
collectWeight(idx + atomStrideX + atomStrideY, mux, muy, 1 - muz);
collectWeight(idx + atomStrideZ, 1 - mux, 1 - muy, muz);
collectWeight(idx + atomStrideX + atomStrideZ, mux, 1 - muy, muz);
collectWeight(idx + atomStrideY + atomStrideZ, 1 - mux, muy, muz);
collectWeight(idx + atomStrideX + atomStrideY + atomStrideZ, mux, muy, muz);
// find dominant atom
let maxWeight = 0.0;
let dominantIdx = -1;
for (const atomIdx in atomWeights) {
if (atomWeights[atomIdx] > maxWeight) {
dominantIdx = atomIdx;
maxWeight = atomWeights[atomIdx];
}
}
if (dominantIdx < 0 || !visibilitySelector.includesAtom(atomLookup[dominantIdx])) {
// this vertex doesn't belong to visible subset and will be skipped
vMap[i] = -1;
continue;
}
}
vMap[i] = newVerCount++;
// color tri-linear interpolation
const dx = (x < xs) ? strideX : 0;
const dy = (y < ys) ? strideY : 0;
const dz = (z < zs) ? strideZ : 0;
const c0 = [0, 0, 0];
const c1 = [0, 0, 0];
const c2 = [0, 0, 0];
const c3 = [0, 0, 0];
idx = colorMap.getDirectIdx(x, y, z);
interp(mux, idx, idx + dx, c0);
interp(mux, idx + dy, idx + dx + dy, c1);
interp(mux, idx + dz, idx + dx + dz, c2);
interp(mux, idx + dy + dz, idx + dx + dy + dz, c3);
const cz0 = [0, 0, 0];
cz0[0] = (1 - muy) * c0[0] + muy * c1[0];
cz0[1] = (1 - muy) * c0[1] + muy * c1[1];
cz0[2] = (1 - muy) * c0[2] + muy * c1[2];
const cz1 = [0, 0, 0];
cz1[0] = (1 - muy) * c2[0] + muy * c3[0];
cz1[1] = (1 - muy) * c2[1] + muy * c3[1];
cz1[2] = (1 - muy) * c2[2] + muy * c3[2];
colors[ind] = (1 - muz) * cz0[0] + muz * cz1[0];
colors[ind + 1] = (1 - muz) * cz0[1] + muz * cz1[1];
colors[ind + 2] = (1 - muz) * cz0[2] + muz * cz1[2];
}
this._colors = colors;
if (visibilitySelector != null) {
// shift visible vertices towards beginning of array
for (i = 0; i < numVerts; ++i) {
const j = vMap[i];
if (j < 0) {
continue;
}
// assert: j <= i
this._position[j * 3] = this._position[i * 3];
this._position[j * 3 + 1] = this._position[i * 3 + 1];
this._position[j * 3 + 2] = this._position[i * 3 + 2];
this._normals[j * 3] = this._normals[i * 3];
this._normals[j * 3 + 1] = this._normals[i * 3 + 1];
this._normals[j * 3 + 2] = this._normals[i * 3 + 2];
this._colors[j * 3] = this._colors[i * 3];
this._colors[j * 3 + 1] = this._colors[i * 3 + 1];
this._colors[j * 3 + 2] = this._colors[i * 3 + 2];
}
// rebuild index list
const numTriangles = this._indices.length / 3;
let newTriCount = 0;
for (i = 0; i < numTriangles; ++i) {
const i0 = vMap[this._indices[3 * i]];
const i1 = vMap[this._indices[3 * i + 1]];
const i2 = vMap[this._indices[3 * i + 2]];
if (i0 >= 0 && i1 >= 0 && i2 >= 0) {
this._indices[3 * newTriCount] = i0;
this._indices[3 * newTriCount + 1] = i1;
this._indices[3 * newTriCount + 2] = i2;
++newTriCount;
}
}
// shrink arrays to data size
this._position = new Float32Array(this._position.buffer.slice(0, newVerCount * 3 * 4));
this._normals = new Float32Array(this._normals.buffer.slice(0, newVerCount * 3 * 4));
this._colors = new Float32Array(this._colors.buffer.slice(0, newVerCount * 3 * 4));
this._indices = new Uint32Array(this._indices.buffer.slice(0, newTriCount * 3 * 4));
}
}