|
@@ -62,6 +62,11 @@ export async function GaussianDensity(ctx: RuntimeContext, unit: Unit, structure
|
|
|
const data = space.create()
|
|
|
const field = Tensor.create(space, data)
|
|
|
|
|
|
+ const idData = space.create()
|
|
|
+ const idField = Tensor.create(space, idData)
|
|
|
+
|
|
|
+ const densData = space.create()
|
|
|
+
|
|
|
const c = Vec3.zero()
|
|
|
|
|
|
const alpha = smoothness
|
|
@@ -74,7 +79,8 @@ export async function GaussianDensity(ctx: RuntimeContext, unit: Unit, structure
|
|
|
const beg = Vec3.zero()
|
|
|
const end = Vec3.zero()
|
|
|
|
|
|
- console.time('density grid')
|
|
|
+ const gridPad = 1 / Math.max(...delta)
|
|
|
+
|
|
|
for (let i = 0; i < elementCount; i++) {
|
|
|
l.element = elements[i]
|
|
|
pos(elements[i], v)
|
|
@@ -85,7 +91,7 @@ export async function GaussianDensity(ctx: RuntimeContext, unit: Unit, structure
|
|
|
const radius = sizeTheme.size(l) + radiusOffset
|
|
|
const rSq = radius * radius
|
|
|
|
|
|
- const r2 = (radiusOffset + radius * 2)
|
|
|
+ const r2 = radiusOffset + radius * 2 + gridPad
|
|
|
const radius2 = Vec3.create(r2, r2, r2)
|
|
|
Vec3.mul(radius2, radius2, delta)
|
|
|
const r2sq = r2 * r2
|
|
@@ -100,7 +106,12 @@ export async function GaussianDensity(ctx: RuntimeContext, unit: Unit, structure
|
|
|
Vec3.div(p, p, delta)
|
|
|
const distSq = Vec3.squaredDistance(p, v)
|
|
|
if (distSq <= r2sq) {
|
|
|
- space.add(data, x, y, z, Math.exp(-alpha * (distSq / rSq)))
|
|
|
+ const dens = Math.exp(-alpha * (distSq / rSq))
|
|
|
+ space.add(data, x, y, z, dens)
|
|
|
+ if (dens > space.get(densData, x, y, z)) {
|
|
|
+ space.set(densData, x, y, z, dens)
|
|
|
+ space.set(idData, x, y, z, i)
|
|
|
+ }
|
|
|
}
|
|
|
}
|
|
|
}
|
|
@@ -110,59 +121,10 @@ export async function GaussianDensity(ctx: RuntimeContext, unit: Unit, structure
|
|
|
await ctx.update({ message: 'filling density grid', current: i, max: elementCount });
|
|
|
}
|
|
|
}
|
|
|
- console.timeEnd('density grid')
|
|
|
-
|
|
|
- const t = Mat4.identity()
|
|
|
- Mat4.fromScaling(t, Vec3.inverse(Vec3.zero(), delta))
|
|
|
- Mat4.setTranslation(t, expandedBox.min)
|
|
|
-
|
|
|
- const { dimensions, get } = space
|
|
|
- const [ xn, yn, zn ] = dimensions
|
|
|
-
|
|
|
- const n = xn * yn * zn
|
|
|
- const idData = space.create()
|
|
|
- const idField = Tensor.create(space, idData)
|
|
|
-
|
|
|
- const lookup3d = unit.lookup3d
|
|
|
|
|
|
- let i = 0
|
|
|
+ const transform = Mat4.identity()
|
|
|
+ Mat4.fromScaling(transform, Vec3.inverse(Vec3.zero(), delta))
|
|
|
+ Mat4.setTranslation(transform, expandedBox.min)
|
|
|
|
|
|
- // TODO use max radius of radiusTheme instead of 4
|
|
|
- const _rSq = 4 + 1 / Math.max(...delta)
|
|
|
- const _minDsq = _rSq * 4
|
|
|
-
|
|
|
- console.time('id grid')
|
|
|
- for (let x = 0; x < xn; ++x) {
|
|
|
- for (let y = 0; y < yn; ++y) {
|
|
|
- for (let z = 0; z < zn; ++z) {
|
|
|
- const dens = get(data, x, y, z)
|
|
|
- let group = -1
|
|
|
- if (dens > 0 && dens < 2) {
|
|
|
- Vec3.set(p, x, y, z)
|
|
|
- Vec3.transformMat4(p, p, t)
|
|
|
- const r = lookup3d.find(p[0], p[1], p[2], _rSq)
|
|
|
- let minDsq = _minDsq
|
|
|
- for (let j = 0, jl = r.count; j < jl; ++j) {
|
|
|
- const dSq = r.squaredDistances[j]
|
|
|
- if (dSq < minDsq) {
|
|
|
- minDsq = dSq
|
|
|
- group = r.indices[j]
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
- space.set(idData, x, y, z, group)
|
|
|
- if (i % 1000000 === 0 && ctx.shouldUpdate) {
|
|
|
- await ctx.update({ message: 'filling id grid', current: i, max: n });
|
|
|
- }
|
|
|
- ++i
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
- console.timeEnd('id grid')
|
|
|
-
|
|
|
- return {
|
|
|
- field,
|
|
|
- idField,
|
|
|
- transform: t
|
|
|
- }
|
|
|
+ return { field, idField, transform }
|
|
|
}
|