diff --git a/scripts/Manifest.toml b/scripts/Manifest.toml
index c1fdc71587300edab541dfdc115a35e364cca3df..0ad0147ff4f18b2616de879323d722c2ffe46422 100644
--- a/scripts/Manifest.toml
+++ b/scripts/Manifest.toml
@@ -1,6 +1,6 @@
# This file is machine-generated - editing it directly is not advised
-julia_version = "1.7.0"
+julia_version = "1.7.1"
manifest_format = "2.0"
[[deps.AbstractFFTs]]
@@ -54,9 +54,9 @@ version = "0.9.11"
[[deps.Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
+git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.16.1+0"
+version = "1.16.1+1"
[[deps.Calculus]]
deps = ["LinearAlgebra"]
@@ -361,9 +361,9 @@ version = "0.21.0+0"
[[deps.Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "74ef6288d071f58033d54fd6708d4bc23a8b8972"
+git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.68.3+1"
+version = "2.68.3+2"
[[deps.Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
@@ -725,9 +725,9 @@ version = "1.10.8"
[[deps.Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
+git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
-version = "1.3.5+0"
+version = "1.3.5+1"
[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
@@ -1289,9 +1289,9 @@ version = "1.6.38+0"
[[deps.libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
-git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
+git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
-version = "1.3.7+0"
+version = "1.3.7+1"
[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index d288ef042da62b86f8a8cc54f9a0976f6fc3ce1b..477c081fb7eea1047bf4128740a7f992bc69384d 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -1,4 +1,6 @@
-using LinearAlgebra: norm, dot, I
+using LinearAlgebra: I, det, dot, norm, normalize
+using SparseArrays: sparse
+using Statistics: mean
using Colors: Gray
# avoid world-age-issues by preloading ColorTypes
@@ -6,13 +8,16 @@ import ColorTypes
import CSV
import DataFrames: DataFrame
import FileIO
+using ForwardDiff: jacobian
using ImageQualityIndexes: assess_psnr, assess_ssim
using OpticalFlowUtils
using WriteVTK: paraview_collection
using Plots
+using StaticArrays: MVector, SArray, SMatrix, SVector, SA
using SemiSmoothNewton
-using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs
+using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs, diam,
+ elmap
using SemiSmoothNewton: project!, project_l2_lagrange!, project_qi_lagrange!,
project_l2_pixel!
using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
@@ -125,6 +130,7 @@ function OptFlowState(mesh;
# tdata will be something like nabla(fw)
T(tdata, u) = tdata * u
+ T(::typeof(adjoint), tdata, v) = tdata' * v
S(u, nablau) = nablau
est = FeFunction(Vest, name = "est")
@@ -446,9 +452,29 @@ end
huber(x, gamma) = abs(x) < gamma ? x^2 / (2 * gamma) : abs(x) - gamma / 2
+# TODO: finish!
+function refine_and_estimate_pd(st::L1L2TVState)
+ # globally refine
+ marked_cells = Set(axes(st.mesh.cells, 2))
+ mesh_new, fs_new = refine(st.mesh, marked_cells;
+ st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
+ st_new = L1L2TVState(mesh_new, st.d, st.m, T, nothing, S,
+ st.alpha1, st.alpha2, st.beta, st.lambda,
+ st.gamma1, st.gamma2,
+ fs_new.est, fs_new.g, fs_new.u, fs_new.p1, fs_new.p2,
+ fs_new.du, fs_new.dp1, fs_new.dp2)
+
+ # compute primal dual error indicators
+ estimate_pd!(st_new)
+
+ # transfer data to old state
+ #
+end
+
# this computes the primal-dual error indicator which is not really useful
-# if not computed on a finer mesh than `u` was solved on
-function estimate!(st::L1L2TVState)
+# if not computed on a finer mesh than `u` was solved on and has problems with
+# alpha1 > 0
+function estimate_pd!(st::L1L2TVState)
function estf(x_; g, u, p1, p2, nablau, w, nablaw, tdata)
alpha1part =
st.alpha1 * huber(norm(st.T(tdata, u) - g), st.gamma1) -
@@ -482,6 +508,95 @@ function estimate!(st::L1L2TVState)
st.est.data .= sqrt.(st.est.data)
end
+# this computes the residual error indicator which has missing theoretical
+# support for alpha1 > 0
+# TODO: finish!
+function estimate_res!(st::L1L2TVState)
+ # FIXME: dirty hack
+ isSnabla = st.S(0, 1) == 1
+ norm2(v) = dot(v, v)
+
+ cellf(hcell; g, u, nablau, p1, p2, tdata) =
+ st.alpha2 *
+ st.T(adjoint, tdata, st.alpha2 * st.T(tdata, u) - g) +
+ st.T(adjoint, tdata, p1)
+ facetf(hfacet, n; g, u, nablau, p1, p2, tdata) =
+ # norm2 is calculated later after both facet contributions are
+ # consolidated
+ n' * (st.beta * st.S(zero(u), nablau) + p2)
+
+ # manual method
+
+ mesh = st.mesh
+ fs = (;st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
+ space = st.est.space
+
+ facetV = Dict{Tuple{Int, Int}, MVector{st.m, Float64}}()
+
+ for cell in cells(mesh)
+ A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
+ view(mesh.vertices, :, view(mesh.cells, :, cell)))
+
+ for f in fs
+ bind!(f, cell)
+ end
+
+ # cell contribution
+
+ hcell = diam(mesh, cell)
+ delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+ intel = abs(det(delmap))
+
+ centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(A, dims = 2))
+ lcentroid = SA[1/3, 1/3]
+
+ opvalues = map(f -> evaluate(f, lcentroid), fs)
+
+ cellres = (isSnabla ? hcell^2 : hcell) *
+ norm2(cellf(hcell; opvalues...)) .* intel
+
+ gdofs = space.dofmap[:, 1, cell]
+ st.est.data[gdofs] .= cellres
+
+ # facet contributions
+
+ cross(v) = SA[-v[2], v[1]]
+ for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
+ fi = mesh.cells[i, cell]
+ fj = mesh.cells[j, cell]
+
+ v = A[:, j] - A[:, i]
+ hfacet = norm(v)
+ normal = -normalize(cross(v))
+ intel = hfacet
+
+ # we don't need to recompute values since only piecewise
+ # constant data is used in `facetf`
+ # we use sqrt here since this value will be squared later
+ facetres = (isSnabla ? sqrt(hfacet) : inv(sqrt(hfacet))) *
+ facetf(hfacet, normal; opvalues...) .* sqrt(intel)
+
+ # average facet contributions
+ v = SVector(facetres)
+ facetV[(fi, fj)] = get(facetV, (fi, fj), zero(v)) + SVector(v)
+ facetV[(fj, fi)] = get(facetV, (fj, fi), zero(v)) + SVector(v)
+ end
+ end
+
+ # add facet contributions to cells
+ for cell in cells(mesh)
+ for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
+ fi = mesh.cells[i, cell]
+ fj = mesh.cells[j, cell]
+
+ gdofs = space.dofmap[:, 1, cell]
+ st.est.data[gdofs] .+= norm2(facetV[(fi, fj)])
+ end
+ end
+
+ st.est.data .= sqrt.(st.est.data)
+end
+
estimate_error(st::L1L2TVState) =
sqrt(sum(x -> x^2, st.est.data) / area(st.mesh))
@@ -596,7 +711,7 @@ function denoise(img; params...)
while true
k += 1
step!(ctx)
- estimate!(ctx)
+ estimate_pd!(ctx)
pvd[k] = save_denoise(ctx, k)
println()
@@ -687,7 +802,7 @@ function denoise_pd(st, img; df=nothing, algorithm, params_...)
domain_factor = 1 / sqrt(area(mesh))
norm_step_ = norm_step(st) * domain_factor
- #estimate!(st)
+ #estimate_pd!(st)
log!(df; k,
norm_step = norm_step_,
@@ -820,7 +935,7 @@ function denoise_approximation(ctx)
k += 1
norm_residual_ = norm_residual(st) * domain_factor
- estimate!(st)
+ estimate_pd!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
@@ -834,7 +949,7 @@ function denoise_approximation(ctx)
if haskey(ctx.params, :tol) && norm_step_ < ctx.params.tol
k += 1
norm_residual_ = norm_residual(st) * domain_factor
- estimate!(st)
+ estimate_pd!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
@@ -858,7 +973,7 @@ function denoise_approximation(ctx)
w = fs.w
norm_residual_ = norm_residual(st) * domain_factor
- estimate!(st)
+ estimate_pd!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
@@ -925,7 +1040,7 @@ function inpaint(img, imgmask; params...)
pvd[0] = save_inpaint(0)
for i in 1:3
step!(ctx)
- estimate!(ctx)
+ estimate_pd!(ctx)
pvd[i] = save_inpaint(i)
println()
end
@@ -1010,7 +1125,7 @@ function optflow(ctx)
k += 1
println()
step!(st)
- #estimate!(st)
+ #estimate_pd!(st)
norm_step_ = norm_step(st) / sqrt(mesh_area)
println("norm_step = $norm_step_")
@@ -1025,7 +1140,7 @@ function optflow(ctx)
display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
- estimate!(st)
+ estimate_res!(st)
pvd[i] = save_step(i)
println("refine ...")