diff --git a/src/graphcut/normalized_cut.jl b/src/graphcut/normalized_cut.jl index 938ad139a..9dabf8d10 100644 --- a/src/graphcut/normalized_cut.jl +++ b/src/graphcut/normalized_cut.jl @@ -117,26 +117,39 @@ end function _recursive_normalized_cut(W, thres, num_cuts) m, n = size(W) - D = Diagonal(vec(sum(W; dims=2))) - - m == 1 && return [1] + (m <= 1) && return ones(Int, m) # trivial + D = Diagonal(vec(sum(W, dims=2))) + + # check that the diagonal is not degenerated as otherwise invDroot errors + dnz = abs.(diag(D)) .>= 1E-16 + if !all(dnz) + # vertices with incident edges summing to almost zero + # are not connected to the rest of the subnetwork, + # put them to separate modules and cut the remaining submatrix + nzlabels = _recursive_normalized_cut(W[dnz, dnz], thres, num_cuts) + nzix = 0 + zix = maximum(nzlabels) + return Int[nz ? nzlabels[nzix += 1] : (zix += 1) for nz in dnz] + end - # get eigenvector corresponding to second smallest eigenvalue + # get eigenvector corresponding to the second smallest generalized eigenvalue: # v = eigs(D-W, D, nev=2, which=SR())[2][:,2] # At least some versions of ARPACK have a bug, this is a workaround invDroot = sqrt.(inv(D)) # equal to Cholesky factorization for diagonal D if n > 12 - λ, Q = eigs(invDroot' * (D - W) * invDroot; nev=12, which=SR()) - ret = real(Q[:, 2]) + _, Q = eigs(invDroot' * (D - W) * invDroot, nev=12, which=SR()) + (size(Q, 2) <= 1) && return collect(1:m) # no 2nd eigenvector + ret = convert(Vector, real(view(Q, :, 2))) else ret = eigen(Matrix(invDroot' * (D - W) * invDroot)).vectors[:, 2] end - v = invDroot * ret + v = real(invDroot * ret) # perform n-cuts with different partitions of v and find best one min_cost = Inf best_thres = -1 - for t in range(minimum(v); stop=maximum(v), length=num_cuts) + vmin, vmax = extrema(v) + for t in range(vmin, stop=vmax, length=num_cuts) cut = v .> t cost = _normalized_cut_cost(cut, W, D) if cost < min_cost diff --git a/src/linalg/LinAlg.jl b/src/linalg/LinAlg.jl index fb5b75e51..3d4cd31d7 100644 --- a/src/linalg/LinAlg.jl +++ b/src/linalg/LinAlg.jl @@ -51,13 +51,8 @@ function eigs(A; kwargs...) schr = partialschur(A; kwargs...) vals, vectors = partialeigen(schr[1]) reved = (kwargs[:which] == LR() || kwargs[:which] == LM()) - k::Int = get(kwargs, :nev, length(vals)) - k = min(k, length(vals)) - perm = collect(1:k) - if vals[1] isa (Real) - perm = sortperm(vals; rev=reved) - perm = perm[1:k] - end + k = min(get(kwargs, :nev, length(vals))::Int, length(vals)) + perm = sortperm(vals, by=real, rev=reved)[1:k] λ = vals[perm] Q = vectors[:, perm] return λ, Q