Skip to content

Mistake in 3D curvature function #61

@Nul1-Ptr

Description

@Nul1-Ptr

Description

in src

function curvature(ϕ::LevelSet, I)
    N = dimension(ϕ)
    if N == 2
        ϕx  = D⁰(ϕ, I, 1)
        ϕy  = D⁰(ϕ, I, 2)
        ϕxx = D2⁰(ϕ, I, 1)
        ϕyy = D2⁰(ϕ, I, 2)
        ϕxy = D2(ϕ, I, (2, 1))
        κ   = (ϕxx * (ϕy)^2 - 2 * ϕy * ϕx * ϕxy + ϕyy * ϕx^2) / (ϕx^2 + ϕy^2)^(3 / 2)
        return κ
    elseif N == 3
        ϕx  = D⁰(ϕ, I, 1)
        ϕy  = D⁰(ϕ, I, 2)
        ϕz  = D⁰(ϕ, I, 3)
        ϕxx = D2⁰(ϕ, I, 1)
        ϕyy = D2⁰(ϕ, I, 2)
        ϕzz = D2⁰(ϕ, I, 3)
        ϕxy = D2(ϕ, I, (2, 1))
        ϕxz = D2(ϕ, I, (3, 1))
        κ   = (ϕxx * (ϕy)^2 - 2 * ϕy * ϕx * ϕxy + ϕyy * ϕx^2 + ϕx^2 * ϕzz - 2 * ϕx * ϕz * ϕxz + ϕz^2 * ϕxx + ϕy^2 * ϕzz - 2 * ϕy * ϕz * ϕyz + ϕz^2 * ϕyy) / (ϕx^2 + ϕy^2)^3 / 2
        return κ
    else
        # generic method
        ∇ϕ = gradient(ϕ, I)
        nrmsq = dot(∇ϕ, ∇ϕ)
        nrmsq == 0.0 && return 0.0= hessian(ϕ, I)
        Δϕ = tr(Hϕ)
        κ = (Δϕ * nrmsq - ∇ϕ' * Hϕ * ∇ϕ) / nrmsq^(3 / 2)
        return κ
    end
end

besides adding $\Phi_{yz}$ in numerator, should the denominator be $(\Phi_x ^2 + \Phi_y ^2 + \Phi_z ^2)^{3/2}$ ?

    elseif N == 3
        ϕx  = D⁰(ϕ, I, 1)
        ϕy  = D⁰(ϕ, I, 2)
        ϕz  = D⁰(ϕ, I, 3)
        ϕxx = D2⁰(ϕ, I, 1)
        ϕyy = D2⁰(ϕ, I, 2)
        ϕzz = D2⁰(ϕ, I, 3)
        ϕxy = D2(ϕ, I, (2, 1))
        ϕxz = D2(ϕ, I, (3, 1))
        ϕyz = D2(ϕ, I, (3, 2))
        κ   = (ϕxx * ϕy^2 + ϕyy * ϕx^2 + ϕxx * ϕz^2 + ϕzz * ϕx^2 + ϕyy * ϕz^2 + ϕzz  * ϕy^2 - 2 * ϕx * ϕz * ϕxz - 2 * ϕy * ϕz * ϕyz  - 2 * ϕy * ϕx * ϕxy) / (ϕx^2 + ϕy^2 + ϕz^2)^(3 / 2)
        return κ

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions