Reconstruct 3-D shape from random Extended Gauss Image

Nanfang Hong, Advanced Study Program

1

A convex 3-D object is one-to-one to Extended Gauss Image (EGI) up to translation

Objective: find height $H$ (distance from origin to each surface)

Define a type ExtendedGaussImage.

In [2]:
using LinearAlgebra
In [3]:
struct ExtendedGaussImage 
    Normal::Array{Float64,2} where T <: Union{Int,Float64}
    Area::Union{Array{T,2},Array{T,1}} where T <: Union{Int,Float64}
    function ExtendedGaussImage((Normal, Area)) 
        if size(Normal)[2] != length(Area)
            error("Number of surface normal does not match number of surface area. ") 
        end
        if sum([norm(Normal[:,i]) for i in 1:length(Area)] - ones(length(Area))) > 1e-10
            error("Normal not unit vector")
        end
        if sum(Normal*Area) > 1e-10
            error("It is not a valid convex polyhedron (Does not satisfy sum of surfaces normal times area equals 0")
        end
        new(Normal,Area) 
    end
end

Generate random extended gauss image

  1. Generate a random matrix N with size(3, faces), then generally rank is 3. Note that to be a polyhedron the faces should at least 4, so the degree of freedom is (faces - 3).
  2. Generate a random vector v2 with length (faces - 3), then can determine v1 such that v = [v1 v2] and N*v = 0.
In [4]:
function randEGI(faces::Int64)
    for i = 1:(10*faces)
        global Normed = zeros(3,1)
        for j = 1:faces  
            N = rand(3)*2 .-1
            N = N ./ norm(N)
            Normed = [Normed N]
        end
        global Normed = Normed[:,2:end]
        
        v2 = rand(faces - 3, 1)
        v1 = inv(Normed[:, 1:3])*(-Normed[:, 4:end]*v2)
        global v = [v1;v2]
        if (sum(v .> zeros(faces,1)) == faces)
            break
        end
    end
    return Normed,v
end
Out[4]:
randEGI (generic function with 1 method)

Example: Cube

In [ ]:
faces = 6
img = ExtendedGaussImage((Float64.([1 -1 0 0 0 0;
    0 0 1 -1 0 0;
    0 0 0 0 1 -1]),Float64.([1,1,1,1,1,1][:,:])))

Example Tetrahedron

In [ ]:
faces = 4
img = ExtendedGaussImage((Float64.([1 0 0 -1;
    0 1 0 -1 ;
    0 0 1 -1 ]),Float64.([1/2,1/2,1/2,3/2][:,:])))

Little's Algorithm (original version)

Little's Algorithm (1985) provides an iterative scheme to reconstruct 3D shape from $EGI = (N, A)$, where $N$ is surface normal of each face and $A$ is area of each face

Let $H$ be distance from origin to each surface, then given $H$ and surface normal $N$, we can uniquely determine a corresponding polyhedron $P = P(H,N)$ by,

  1. Map the $n$ planes given by $H$ into $M$, a set of $n$ points in $R^3$, using the dual transform. Dual transform is defined as, given a plane $$Ax + By + Cz + 1 = 0$$, then corresponding dual point is $$(A, B, C)$$

  2. Compute convex hull of $M$, which is also called the dual of $P$ . Then any face of the dual of $P$ corresponds a vertex of $P$, any two points incident on an edge in the dual of $P$ correspond to a pair of faces of $P$ which share an edge. Then we know $P$ by its vertices and faces.

Little's Algorithm is:

1. Initial $H$ as

$$[1,1,...,1]$$

2. Compute $P$

3. Compute $V(P)$ by area of faces $S = [S_1, S_2, ... ,S_n]$ using

$$V(P) = \frac{1}{3}H^TS$$

gradient of V(P(H,N)) with respect to H is constant $$\nabla V(P) = \frac{1}{3}S$$

updating H by multiple of step $\Delta$,

$$\Delta = \langle\,A,\nabla V(P)\rangle \nabla V(P) - A $$

which means just update $H$ proportional to $A$. Before update, scale $H$ by $V(P)^\frac{1}{3}=(\frac{1}{3}H^TS)^\frac{1}{3}$ to keep its volume unity

$$H \leftarrow(\frac{1}{3}H^TS)^\frac{1}{3}H - k\Delta$$

4. repeat 2 and 3 until $|(H_{t+1} - H_{t})A|$ is less than $\epsilon$

In [5]:
using Polyhedra, QHull, Makie
In [6]:
function Surface(p::QHull.Polyhedron)
    S = Float64[]
    for pi in eachindex(halfspaces(p))
        global n = get(p, pi).a
        n = n/norm(n)
        get(p, pi).β
        global pts = incidentpoints(p, pi)
        global b = pts[2] - pts[1]
        b = b/norm(b)
        global a = cross(n, b) 
        norm(a)
        append!(S, Float64(Polyhedra.volume(polyhedron(vrep([[a b]'*pt for pt in pts]), QHull.Library()))))
    end
    return S
end
Out[6]:
Surface (generic function with 1 method)
In [7]:
function Volume(H::Array{Float64,1}, S::Array{Float64,1}) 
    return 1/3*(H'*S)
end

function Volume(p::QHull.Polyhedron) 
    faces = length(hrep(p).b)
    H = MixedMatHRep(hrep(p)).b ./[norm(MixedMatHRep(hrep(p)).A[i,:]) for i in 1:faces]
    S = Float64[]
    for pi in eachindex(halfspaces(p))
        global n = get(p, pi).a
        n = n/norm(n)
        get(p, pi).β
        global pts = incidentpoints(p, pi)
        global b = pts[2] - pts[1]
        b = b/norm(b)
        global a = cross(n, b) 
        norm(a)
        append!(S, Float64(Polyhedra.volume(polyhedron(vrep([[a b]'*pt for pt in pts]), QHull.Library()))))
    end
    return 1/3*(H'*S)
end
Out[7]:
Volume (generic function with 2 methods)
In [8]:
function Polyhedron(H,N)
    h = hrep(N', [norm(N[:,i]) for i in 1:length(H)] .* H)
    p = polyhedron(h, QHull.Library())
    return p
end
Out[8]:
Polyhedron (generic function with 1 method)
In [9]:
function Normal(p::QHull.Polyhedron) 
    faces = length(hrep(p).b)
    return [MixedMatHRep(hrep(p)).A[i,:]./norm(MixedMatHRep(hrep(p)).A[i,:]) for i in 1:faces]
end
Out[9]:
Normal (generic function with 1 method)
In [10]:
function randPolyhedron(faces)
    egi = randEGI(faces)
    img = ExtendedGaussImage(egi)
    for i in 1:(10*faces)
        H = rand(faces)
        global p = Polyhedron(H, img.Normal)
        if length(Surface(p)) == faces
            break
        end
    end
    return p, ExtendedGaussImage((img.Normal, Surface(p)[:,:]))
end
Out[10]:
randPolyhedron (generic function with 1 method)

Given a random valid EGI

In [11]:
faces = 6#if faces = 4, it reconstruct immediately
p, img = randPolyhedron(faces)
Out[11]:
(HalfSpace([-5.67972, -2.81342, 2.77206], 1.0) ∩ HalfSpace([-1.36039, 1.36022, 0.715521], 1.0) ∩ HalfSpace([-0.80731, 0.0541164, -0.847875], 1.0) ∩ HalfSpace([0.83325, 0.254902, 0.873617], 1.0) ∩ HalfSpace([0.0901262, -2.06049, -0.509373], 1.0) ∩ HalfSpace([0.949222, 0.733648, 0.645261], 1.0) : convexhull([-0.6361, 0.387788, -0.549001], [0.510115, 2.05235, -1.53414], [-0.416483, -0.305192, -0.802341], [3.48488, 0.766853, -4.44863], [0.108731, 0.35007, 0.938817], [0.134997, 0.396668, 0.900169], [0.54938, -0.662514, 0.813979], [1.43146, -0.396782, -0.104877]), ExtendedGaussImage([-0.821005 -0.662789 … 0.0424237 0.696824; -0.40668 0.66271 … -0.969903 0.538571; 0.400702 0.348606 … -0.239769 0.473686], [1.65763; 1.87756; … ; 6.06886; 9.02011]))
In [12]:
function RenderPolyhedron(p::QHull.Polyhedron)
    m = Polyhedra.Mesh(p)
    scene = mesh(m, color=:white)
end
Out[12]:
RenderPolyhedron (generic function with 1 method)
In [13]:
RenderPolyhedron(p)
Out[13]:

We now use Little's algorithm to recontruct $q$ from the EGI. We want $p = q$

In [14]:
H = ones(faces)
q = Polyhedron(H, img.Normal)
RenderPolyhedron(q)
Out[14]:

$p$ and $q$ have same surface normal

In [17]:
Normal(q)
Out[17]:
6-element Array{Array{Float64,1},1}:
 [-0.821005, -0.40668, 0.400702]  
 [-0.662789, 0.66271, 0.348606]   
 [-0.688835, 0.0461746, -0.723446]
 [0.675303, 0.206584, 0.708018]   
 [0.0424237, -0.969903, -0.239769]
 [0.696824, 0.538571, 0.473686]   
In [18]:
Normal(p)
Out[18]:
6-element Array{Array{Float64,1},1}:
 [-0.821005, -0.40668, 0.400702]  
 [-0.662789, 0.66271, 0.348606]   
 [-0.688835, 0.0461746, -0.723446]
 [0.675303, 0.206584, 0.708018]   
 [0.0424237, -0.969903, -0.239769]
 [0.696824, 0.538571, 0.473686]   

Little's Algorithm (slightly modified)

Use KL divergence of proposed and true area distributions as error measure

Modify iterative formula


Given $EGI(N,A)$ of a convex polyhedron $p$ , want to construct polyhedron $q$ such that surface area distribution $P_A$ and $P_S$ are the same (KL divergence 0)

($N$ is surface normal of $p$ and $q$, $A$ is surface area of $p$, $S$ is surface area of $q$. )

Initialize the $H$ (distance from origin to every surface) of $q$ all be 1, adjust $H$ such that KL divergence converges.

1. Initial $H$ as $[1,1,...,1]$, $P_A \propto A$

2. Compute $q$

3. Scaling $H$ such that $V(q) = 1$, which means $H \leftarrow \frac{H}{V(q)^\frac{1}{3}}$

4. Compute $q$ again. Compute $S$

5. $\nabla V(q) = \frac{1}{3}S$

6. $\Delta = \nabla V(q) - \frac{A}{\mid A \mid} {\mid V(q) \mid} $

7. $H \leftarrow H + k\Delta$

8. Scaling $H$ again, $H \leftarrow \frac{H}{V(q)^\frac{1}{3}}$

9. Compute $q$ again. Compute $S$. $P_S \propto S$

10. Repeat 3-9 until $D(P_S \mid \mid P_A) < \epsilon$

In [ ]:
Surface(q)
In [15]:
Base.length(q::QHull.Polyhedron) = 1
Base.iterate(q::QHull.Polyhedron) = (q, nothing)
function Little(k, ϵ, img::ExtendedGaussImage)
    faces = length(img.Area)
    H = ones(faces)
    p1 = img.Area / sum(img.Area)
    KLArray = []
    KL = 0.05
    qArray = QHull.Polyhedron[]
    q = Polyhedron(H, img.Normal)
    while KL >= ϵ
        H = H / (Volume(q))^(1/3)
        q = Polyhedron(H, img.Normal)
         = 1/3*Surface(q)
        Δ =  - (img.Area * norm())/ norm(img.Area)
        H = vec(H + k .* Δ)
        H = H / (Volume(q))^(1/3)
        q = Polyhedron(H, img.Normal)
        p2 = Surface(q) / sum(Surface(q))
        KL = sum([p2[i]*log(p2[i]/p1[i]) for i in 1:faces])
        @show p2 ./ p1
        @show KL 
        append!(KLArray, KL)
        append!(qArray, q)
    end
    return qArray, KLArray
end
Out[15]:
Little (generic function with 1 method)
In [16]:
k = 0.5
ϵ = 1e-11
qArray, KL = Little(k, ϵ, img)
p2 ./ p1 = [0.987372; 0.957916; 1.02358; 1.46553; 0.97987; 0.967926]
KL = 0.003188198216623906
p2 ./ p1 = [1.01082; 0.964599; 1.0121; 1.18764; 0.987789; 0.988877]
KL = 0.0006018150551309266
p2 ./ p1 = [1.01265; 0.974821; 1.00682; 1.08633; 0.992406; 0.995754]
KL = 0.00015202446296400563
p2 ./ p1 = [1.00971; 0.983443; 1.00398; 1.04349; 0.995314; 0.99824]
KL = 4.590537225336415e-5
p2 ./ p1 = [1.00657; 0.989534; 1.00236; 1.02336; 0.997134; 0.999209]
KL = 1.5208404262194745e-5
p2 ./ p1 = [1.00421; 0.993523; 1.00141; 1.0131; 0.998257; 0.999616]
KL = 5.268648872459891e-6
p2 ./ p1 = [1.00262; 0.996038; 1.00084; 1.00755; 0.998945; 0.999801]
KL = 1.8650043537053675e-6
p2 ./ p1 = [1.00161; 0.997593; 1.00051; 1.00443; 0.999363; 0.999891]
KL = 6.6735125859043e-7
p2 ./ p1 = [1.00098; 0.998543; 1.0003; 1.00263; 0.999615; 0.999939]
KL = 2.40152082429252e-7
p2 ./ p1 = [1.00059; 0.99912; 1.00018; 1.00157; 0.999768; 0.999964]
KL = 8.668663797134696e-8
p2 ./ p1 = [1.00036; 0.999469; 1.00011; 1.00094; 0.99986; 0.999979]
KL = 3.134447860114585e-8
p2 ./ p1 = [1.00022; 0.99968; 1.00007; 1.00056; 0.999916; 0.999988]
KL = 1.1344640291814223e-8
p2 ./ p1 = [1.00013; 0.999807; 1.00004; 1.00034; 0.999949; 0.999993]
KL = 4.108288273178381e-9
p2 ./ p1 = [1.00008; 0.999884; 1.00002; 1.0002; 0.99997; 0.999996]
KL = 1.4882272288808716e-9
p2 ./ p1 = [1.00005; 0.99993; 1.00001; 1.00012; 0.999982; 0.999997]
KL = 5.392085637500465e-10
p2 ./ p1 = [1.00003; 0.999958; 1.00001; 1.00007; 0.999989; 0.999998]
KL = 1.9538448454964382e-10
p2 ./ p1 = [1.00002; 0.999975; 1.00001; 1.00004; 0.999993; 0.999999]
KL = 7.080272191572292e-11
p2 ./ p1 = [1.00001; 0.999985; 1.0; 1.00003; 0.999996; 0.999999]
KL = 2.5658191261448877e-11
p2 ./ p1 = [1.00001; 0.999991; 1.0; 1.00002; 0.999998; 1.0]
KL = 9.298322946510732e-12
Out[16]:
(QHull.Polyhedron[HalfSpace([-2.26988, -1.12437, 1.10785], 1.0) ∩ HalfSpace([-1.79709, 1.79688, 0.945215], 1.0) ∩ HalfSpace([-1.7014, 0.11405, -1.78689], 1.0) ∩ HalfSpace([1.55761, 0.476493, 1.63306], 1.0) ∩ HalfSpace([0.115461, -2.63971, -0.652562], 1.0) ∩ HalfSpace([2.09215, 1.61701, 1.4222], 1.0) : convexhull([-0.379201, -0.340947, -0.220333], [1.82248, 0.264046, -2.27806], [-0.512355, 0.0792071, -0.0667324], [0.121175, 0.999238, -0.611232], [0.134371, -0.530944, 0.639102], [0.743313, -0.347496, 0.00477013], [-0.12256, 0.0594681, 0.711892], [-0.0710295, 0.150892, 0.636066]), HalfSpace([-2.27377, -1.1263, 1.10974], 1.0) ∩ HalfSpace([-1.8145, 1.81428, 0.954369], 1.0) ∩ HalfSpace([-1.6536, 0.110846, -1.73668], 1.0) ∩ HalfSpace([1.49342, 0.456858, 1.56577], 1.0) ∩ HalfSpace([0.116399, -2.66114, -0.657859], 1.0) ∩ HalfSpace([2.15482, 1.66545, 1.4648], 1.0) : convexhull([-0.385711, -0.335796, -0.229985], [1.80763, 0.266905, -2.27992], [-0.516361, 0.0764591, -0.0792722], [0.117658, 0.9972, -0.624192], [0.144007, -0.531766, 0.656469], [0.686802, -0.368245, 0.0910422], [-0.11176, 0.0559703, 0.728929], [-0.0847227, 0.103938, 0.689145]), HalfSpace([-2.26576, -1.12233, 1.10583], 1.0) ∩ HalfSpace([-1.83013, 1.82992, 0.962593], 1.0) ∩ HalfSpace([-1.63135, 0.109354, -1.71332], 1.0) ∩ HalfSpace([1.46838, 0.449196, 1.53951], 1.0) ∩ HalfSpace([0.117187, -2.67916, -0.662313], 1.0) ∩ HalfSpace([2.17546, 1.6814, 1.47883], 1.0) : convexhull([-0.390378, -0.332679, -0.233196], [1.80302, 0.27004, -2.28319], [-0.518759, 0.0724173, -0.0851005], [0.118563, 0.997954, -0.632858], [0.146344, -0.53124, 0.664979], [0.664491, -0.375144, 0.125228], [-0.107283, 0.0515804, 0.736833], [-0.0881683, 0.0854926, 0.708706]), HalfSpace([-2.2579, -1.11843, 1.10199], 1.0) ∩ HalfSpace([-1.84193, 1.84171, 0.968799], 1.0) ∩ HalfSpace([-1.61986, 0.108584, -1.70125], 1.0) ∩ HalfSpace([1.45707, 0.445736, 1.52765], 1.0) ∩ HalfSpace([0.117776, -2.69263, -0.665644], 1.0) ∩ HalfSpace([2.18285, 1.68711, 1.48385], 1.0) : convexhull([-0.393441, -0.330674, -0.234291], [1.80118, 0.272379, -2.28543], [-0.520145, 0.0691253, -0.0881315], [0.120061, 0.998849, -0.638368], [0.146644, -0.53048, 0.669513], [0.654264, -0.377556, 0.140728], [-0.105142, 0.0481098, 0.740845], [-0.0887562, 0.0771813, 0.716734]), HalfSpace([-2.25217, -1.1156, 1.0992], 1.0) ∩ HalfSpace([-1.84995, 1.84973, 0.973015], 1.0) ∩ HalfSpace([-1.61353, 0.10816, -1.6946], 1.0) ∩ HalfSpace([1.45143, 0.444011, 1.52174], 1.0) ∩ HalfSpace([0.118175, -2.70175, -0.667899], 1.0) ∩ HalfSpace([2.18565, 1.68927, 1.48576], 1.0) : convexhull([-0.395379, -0.329411, -0.234671], [1.80034, 0.273944, -2.28683], [-0.520958, 0.0668416, -0.0898085], [0.121283, 0.999522, -0.641794], [0.146454, -0.529863, 0.672057], [0.649113, -0.378434, 0.148441], [-0.104027, 0.0457285, 0.74302], [-0.0886181, 0.0730665, 0.720346]), HalfSpace([-2.24839, -1.11373, 1.09735], 1.0) ∩ HalfSpace([-1.85511, 1.85489, 0.975729], 1.0) ∩ HalfSpace([-1.60991, 0.107917, -1.6908], 1.0) ∩ HalfSpace([1.44842, 0.44309, 1.51858], 1.0) ∩ HalfSpace([0.118432, -2.70762, -0.669349], 1.0) ∩ HalfSpace([2.18677, 1.69015, 1.48652], 1.0) : convexhull([-0.396579, -0.328629, -0.234805], [1.79992, 0.27494, -2.28769], [-0.521441, 0.0653608, -0.0907698], [0.12213, 0.999972, -0.643898], [0.146214, -0.529436, 0.673528], [0.646346, -0.378767, 0.152543], [-0.103413, 0.0441921, 0.744249], [-0.0883673, 0.0708861, 0.722109]), HalfSpace([-2.24599, -1.11254, 1.09619], 1.0) ∩ HalfSpace([-1.85834, 1.85811, 0.977427], 1.0) ∩ HalfSpace([-1.6078, 0.107775, -1.68858], 1.0) ∩ HalfSpace([1.44673, 0.442575, 1.51682], 1.0) ∩ HalfSpace([0.118592, -2.71128, -0.670254], 1.0) ∩ HalfSpace([2.18726, 1.69052, 1.48685], 1.0) : convexhull([-0.397313, -0.32815, -0.234853], [1.7997, 0.27556, -2.28822], [-0.521729, 0.064432, -0.0913322], [0.122679, 1.00026, -0.64518], [0.146025, -0.529159, 0.674394], [0.644795, -0.378901, 0.154828], [-0.103063, 0.043231, 0.744962], [-0.0881582, 0.069675, 0.723029]), HalfSpace([-2.24451, -1.11181, 1.09546], 1.0) ∩ HalfSpace([-1.86032, 1.8601, 0.978471], 1.0) ∩ HalfSpace([-1.60654, 0.107691, -1.68727], 1.0) ∩ HalfSpace([1.44576, 0.442279, 1.5158], 1.0) ∩ HalfSpace([0.11869, -2.71353, -0.670811], 1.0) ∩ HalfSpace([2.18749, 1.69069, 1.48701], 1.0) : convexhull([-0.397759, -0.327859, -0.23487], [1.79957, 0.27594, -2.28854], [-0.521902, 0.0638598, -0.0916652], [0.123022, 1.00044, -0.645957], [0.145896, -0.528986, 0.674908], [0.6439, -0.378958, 0.156139], [-0.102859, 0.0426397, 0.745381], [-0.0880121, 0.0689811, 0.723534]), HalfSpace([-2.24361, -1.11136, 1.09502], 1.0) ∩ HalfSpace([-1.86153, 1.86131, 0.979108], 1.0) ∩ HalfSpace([-1.6058, 0.107641, -1.68648], 1.0) ∩ HalfSpace([1.4452, 0.442105, 1.51521], 1.0) ∩ HalfSpace([0.11875, -2.7149, -0.67115], 1.0) ∩ HalfSpace([2.1876, 1.69078, 1.48708], 1.0) : convexhull([-0.39803, -0.327683, -0.234877], [1.7995, 0.276171, -2.28873], [-0.522006, 0.0635106, -0.0918638], [0.123234, 1.00055, -0.646427], [0.145813, -0.528879, 0.675215], [0.643376, -0.378985, 0.156907], [-0.102739, 0.0422792, 0.745631], [-0.087917, 0.0685756, 0.723821]), HalfSpace([-2.24306, -1.11108, 1.09475], 1.0) ∩ HalfSpace([-1.86227, 1.86204, 0.979494], 1.0) ∩ HalfSpace([-1.60535, 0.107612, -1.68602], 1.0) ∩ HalfSpace([1.44486, 0.442002, 1.51486], 1.0) ∩ HalfSpace([0.118787, -2.71574, -0.671356], 1.0) ∩ HalfSpace([2.18766, 1.69083, 1.48712], 1.0) : convexhull([-0.398193, -0.327577, -0.234879], [1.79946, 0.276311, -2.28885], [-0.522068, 0.0632988, -0.0919827], [0.123363, 1.00061, -0.64671], [0.145761, -0.528814, 0.675399], [0.643065, -0.378997, 0.15736], [-0.102667, 0.0420606, 0.74578], [-0.0878573, 0.0683359, 0.723987]), HalfSpace([-2.24272, -1.11092, 1.09459], 1.0) ∩ HalfSpace([-1.86271, 1.86249, 0.979728], 1.0) ∩ HalfSpace([-1.60509, 0.107594, -1.68574], 1.0) ∩ HalfSpace([1.44466, 0.441941, 1.51465], 1.0) ∩ HalfSpace([0.118809, -2.71624, -0.671481], 1.0) ∩ HalfSpace([2.18769, 1.69085, 1.48714], 1.0) : convexhull([-0.398292, -0.327513, -0.23488], [1.79944, 0.276395, -2.28892], [-0.522106, 0.0631707, -0.092054], [0.123441, 1.00065, -0.646881], [0.14573, -0.528774, 0.675509], [0.642879, -0.379004, 0.157631], [-0.102624, 0.0419284, 0.745869], [-0.0878205, 0.068193, 0.724086]), HalfSpace([-2.24252, -1.11082, 1.09449], 1.0) ∩ HalfSpace([-1.86298, 1.86276, 0.979869], 1.0) ∩ HalfSpace([-1.60492, 0.107583, -1.68557], 1.0) ∩ HalfSpace([1.44454, 0.441904, 1.51452], 1.0) ∩ HalfSpace([0.118822, -2.71654, -0.671556], 1.0) ∩ HalfSpace([2.18771, 1.69087, 1.48716], 1.0) : convexhull([-0.398351, -0.327474, -0.234881], [1.79942, 0.276446, -2.28896], [-0.522129, 0.0630934, -0.0920968], [0.123489, 1.00068, -0.646984], [0.14571, -0.52875, 0.675576], [0.642768, -0.379008, 0.157793], [-0.102599, 0.0418487, 0.745923], [-0.087798, 0.0681076, 0.724144]), HalfSpace([-2.2424, -1.11076, 1.09443], 1.0) ∩ HalfSpace([-1.86314, 1.86292, 0.979954], 1.0) ∩ HalfSpace([-1.60483, 0.107576, -1.68546], 1.0) ∩ HalfSpace([1.44447, 0.441882, 1.51444], 1.0) ∩ HalfSpace([0.11883, -2.71673, -0.671601], 1.0) ∩ HalfSpace([2.18772, 1.69087, 1.48716], 1.0) : convexhull([-0.398387, -0.327451, -0.234881], [1.79942, 0.276477, -2.28899], [-0.522142, 0.0630468, -0.0921226], [0.123517, 1.00069, -0.647046], [0.145698, -0.528736, 0.675616], [0.642702, -0.37901, 0.15789], [-0.102583, 0.0418006, 0.745955], [-0.0877843, 0.0680563, 0.724179]), HalfSpace([-2.24233, -1.11072, 1.0944], 1.0) ∩ HalfSpace([-1.86324, 1.86302, 0.980005], 1.0) ∩ HalfSpace([-1.60477, 0.107572, -1.6854], 1.0) ∩ HalfSpace([1.44442, 0.441869, 1.5144], 1.0) ∩ HalfSpace([0.118835, -2.71684, -0.671628], 1.0) ∩ HalfSpace([2.18772, 1.69088, 1.48717], 1.0) : convexhull([-0.398408, -0.327437, -0.234881], [1.79941, 0.276495, -2.289], [-0.52215, 0.0630187, -0.0921381], [0.123534, 1.0007, -0.647083], [0.145691, -0.528727, 0.67564], [0.642662, -0.379011, 0.157948], [-0.102574, 0.0417716, 0.745974], [-0.0877761, 0.0680255, 0.7242]), HalfSpace([-2.24228, -1.1107, 1.09438], 1.0) ∩ HalfSpace([-1.8633, 1.86307, 0.980036], 1.0) ∩ HalfSpace([-1.60473, 0.10757, -1.68537], 1.0) ∩ HalfSpace([1.4444, 0.441861, 1.51437], 1.0) ∩ HalfSpace([0.118838, -2.7169, -0.671645], 1.0) ∩ HalfSpace([2.18773, 1.69088, 1.48717], 1.0) : convexhull([-0.398421, -0.327428, -0.234881], [1.79941, 0.276506, -2.28901], [-0.522155, 0.0630018, -0.0921474], [0.123545, 1.0007, -0.647106], [0.145687, -0.528722, 0.675654], [0.642637, -0.379012, 0.157984], [-0.102569, 0.0417542, 0.745986], [-0.0877711, 0.068007, 0.724212]), HalfSpace([-2.24226, -1.11069, 1.09436], 1.0) ∩ HalfSpace([-1.86333, 1.86311, 0.980055], 1.0) ∩ HalfSpace([-1.60471, 0.107569, -1.68534], 1.0) ∩ HalfSpace([1.44438, 0.441856, 1.51436], 1.0) ∩ HalfSpace([0.11884, -2.71694, -0.671655], 1.0) ∩ HalfSpace([2.18773, 1.69088, 1.48717], 1.0) : convexhull([-0.398429, -0.327423, -0.234881], [1.79941, 0.276513, -2.28902], [-0.522158, 0.0629916, -0.092153], [0.123551, 1.00071, -0.647119], [0.145684, -0.528719, 0.675663], [0.642623, -0.379012, 0.158005], [-0.102565, 0.0417437, 0.745993], [-0.0877681, 0.0679958, 0.72422]), HalfSpace([-2.24224, -1.11068, 1.09435], 1.0) ∩ HalfSpace([-1.86335, 1.86313, 0.980066], 1.0) ∩ HalfSpace([-1.6047, 0.107568, -1.68533], 1.0) ∩ HalfSpace([1.44437, 0.441853, 1.51435], 1.0) ∩ HalfSpace([0.118841, -2.71697, -0.671661], 1.0) ∩ HalfSpace([2.18773, 1.69089, 1.48717], 1.0) : convexhull([-0.398434, -0.32742, -0.234881], [1.7994, 0.276517, -2.28902], [-0.52216, 0.0629855, -0.0921564], [0.123555, 1.00071, -0.647127], [0.145683, -0.528717, 0.675668], [0.642614, -0.379013, 0.158017], [-0.102563, 0.0417373, 0.745997], [-0.0877663, 0.0679891, 0.724224]), HalfSpace([-2.24223, -1.11068, 1.09435], 1.0) ∩ HalfSpace([-1.86337, 1.86314, 0.980073], 1.0) ∩ HalfSpace([-1.60469, 0.107567, -1.68532], 1.0) ∩ HalfSpace([1.44437, 0.441851, 1.51434], 1.0) ∩ HalfSpace([0.118841, -2.71698, -0.671664], 1.0) ∩ HalfSpace([2.18773, 1.69089, 1.48717], 1.0) : convexhull([-0.398437, -0.327418, -0.234881], [1.7994, 0.27652, -2.28902], [-0.522161, 0.0629818, -0.0921585], [0.123557, 1.00071, -0.647132], [0.145682, -0.528716, 0.675671], [0.642609, -0.379013, 0.158025], [-0.102562, 0.0417335, 0.746], [-0.0877652, 0.0679851, 0.724227]), HalfSpace([-2.24223, -1.11067, 1.09435], 1.0) ∩ HalfSpace([-1.86337, 1.86315, 0.980077], 1.0) ∩ HalfSpace([-1.60469, 0.107567, -1.68532], 1.0) ∩ HalfSpace([1.44436, 0.44185, 1.51434], 1.0) ∩ HalfSpace([0.118842, -2.71699, -0.671666], 1.0) ∩ HalfSpace([2.18773, 1.69089, 1.48717], 1.0) : convexhull([-0.398438, -0.327417, -0.234881], [1.7994, 0.276521, -2.28903], [-0.522162, 0.0629796, -0.0921597], [0.123558, 1.00071, -0.647135], [0.145681, -0.528715, 0.675673], [0.642606, -0.379013, 0.15803], [-0.102561, 0.0417312, 0.746001], [-0.0877646, 0.0679827, 0.724229])], Any[0.0031882, 0.000601815, 0.000152024, 4.59054e-5, 1.52084e-5, 5.26865e-6, 1.865e-6, 6.67351e-7, 2.40152e-7, 8.66866e-8, 3.13445e-8, 1.13446e-8, 4.10829e-9, 1.48823e-9, 5.39209e-10, 1.95384e-10, 7.08027e-11, 2.56582e-11, 9.29832e-12])
In [19]:
typeof(qArray)
Out[19]:
Array{QHull.Polyhedron,1}
In [ ]:
m = Polyhedra.Mesh(qArray[end])
scene = mesh(m, color=:white)
In [ ]:
using Plots
Plots.plot(1:length(KL), KL, label = "KL Divergence")

5

6

Animation of convergence by Makie.jl

2

Thought

In [ ]:
using Makie

for i in 1:length(qArray)
    @show i
    scene = Scene()
    m = Polyhedra.Mesh(qArray[i])
    mesh!(scene, m, color=:white)
    Makie.save(string(i)*".png", scene)
    sleep(1/24)
end
In [1]:
qArray
UndefVarError: qArray not defined

Stacktrace:
 [1] top-level scope at In[1]:1

Try to use JuMP.jl but fail

$$\min_{H} HA $$

Subject to$$V(q) = 1$$

In [ ]:
using JuMP
In [ ]:
m = Model()
@variable(m, H[1:faces] >= 0)
@constraint(m,  Polyhedra.volume(polyhedron(hrep(img.Normal', [norm(img.Normal[:,i]) for i in 1:faces] .* H), QHull.Library())) == 1)
@objective(m, Min, H'*img.Area)
solve(m)