diff --git a/scripts/Manifest.toml b/scripts/Manifest.toml
new file mode 100644
index 0000000000000000000000000000000000000000..d16b6d678d3dd9219c4482adc26b980b8c059ab4
--- /dev/null
+++ b/scripts/Manifest.toml
@@ -0,0 +1,429 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[CSV]]
+deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
+git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
+uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+version = "0.8.5"
+
+[[ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.3.0"
+
+[[CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.0"
+
+[[ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.11.0"
+
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.8"
+
+[[CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "79b9563ef3f2cc5fc6d3046a5ee1a57c9de52495"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.33.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[Crayons]]
+git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
+uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+version = "4.0.4"
+
+[[DataAPI]]
+git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.7.0"
+
+[[DataFrames]]
+deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
+git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
+uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+version = "1.2.2"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.10"
+
+[[DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[DiffResults]]
+deps = ["StaticArrays"]
+git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.0.3"
+
+[[DiffRules]]
+deps = ["NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.3.0"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.5"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[FileIO]]
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
+uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+version = "1.10.1"
+
+[[FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
+git-tree-sha1 = "8c8eac2af06ce35973c3eadb4ab3243076a408e7"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.12.1"
+
+[[FixedPointNumbers]]
+deps = ["Statistics"]
+git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.8.4"
+
+[[Formatting]]
+deps = ["Printf"]
+git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
+uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
+version = "0.4.2"
+
+[[ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.19"
+
+[[Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[InvertedIndices]]
+deps = ["Test"]
+git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
+uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
+version = "1.0.0"
+
+[[IrrationalConstants]]
+git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.0"
+
+[[IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
+[[JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.3.0"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.16.1+1"
+
+[[LightXML]]
+deps = ["Libdl", "XML2_jll"]
+git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+version = "0.9.0"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[LogExpFunctions]]
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.0"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.7"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "1.0.0"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[NaNMath]]
+git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.5"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "1.1.2"
+
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[PooledArrays]]
+deps = ["DataAPI", "Future"]
+git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
+uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
+version = "1.2.1"
+
+[[Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.2.2"
+
+[[PrettyTables]]
+deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
+git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
+uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+version = "1.1.0"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Reexport]]
+git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.1.0"
+
+[[Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.1.3"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[SemiSmoothNewton]]
+deps = ["CSV", "ColorTypes", "Colors", "DataFrames", "FileIO", "ForwardDiff", "LinearAlgebra", "SparseArrays", "StaticArrays", "Statistics", "WriteVTK"]
+path = ".."
+uuid = "a35a6534-5ddf-4107-9885-a23102237ad0"
+version = "0.1.0"
+
+[[SentinelArrays]]
+deps = ["Dates", "Random"]
+git-tree-sha1 = "a3a337914a035b2d59c9cbe7f1a38aaba1265b02"
+uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
+version = "1.3.6"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.0.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[SpecialFunctions]]
+deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
+git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "1.6.1"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.2.12"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
+git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.5.0"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TranscodingStreams]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.9.6"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[WriteVTK]]
+deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams"]
+git-tree-sha1 = "cc6b182732e3e00c3942f4e84fe88f0ec46e89a7"
+uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
+version = "1.10.1"
+
+[[XML2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
+uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
+version = "2.9.12+0"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/scripts/Project.toml b/scripts/Project.toml
new file mode 100644
index 0000000000000000000000000000000000000000..c15713246bd115374fdcdee2408baf90f6c54163
--- /dev/null
+++ b/scripts/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+SemiSmoothNewton = "a35a6534-5ddf-4107-9885-a23102237ad0"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
diff --git a/src/run.jl b/scripts/run.jl
similarity index 89%
rename from src/run.jl
rename to scripts/run.jl
index 2478e3e345d91d5bc2dbc57abf35b522862dbdfb..cfc70b66eb5f91bb5a19e1048a3813d2d8e1fe39 100644
--- a/src/run.jl
+++ b/scripts/run.jl
@@ -1,12 +1,30 @@
-export myrun, denoise, denoise_pd, inpaint, optflow, solve_primal!, estimate!, loadimg, saveimg
+using LinearAlgebra: norm, dot
-using LinearAlgebra: norm
using Colors: Gray
-
# avoid world-age-issues by preloading ColorTypes
import ColorTypes
+import DataFrames: DataFrame
import FileIO
+using WriteVTK: paraview_collection
+
+using SemiSmoothNewton
+using SemiSmoothNewton: project_img!, project!
+using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
+
+include("util.jl")
+isdefined(Main, :Revise) && Revise.track(joinpath(@__DIR__, "util.jl"))
+
+toimg(arr) = Gray.(clamp.(arr, 0., 1.))
+loadimg(x) = reverse(transpose(Float64.(FileIO.load(x))); dims = 2)
+saveimg(io, x) = FileIO.save(io, transpose(reverse(toimg(x); dims = 2)))
+"""
+ logfilter(dt; a=20)
+
+Reduce dataset such that the resulting dataset would have approximately
+equidistant datapoints on a log-scale. `a` controls density.
+"""
+logfilter(dt; a=20) = filter(:k => k -> round(a*log(k+1)) - round(a*log(k)) > 0, dt)
struct L1L2TVContext{M, Ttype, Stype}
name::String
@@ -72,6 +90,24 @@ function L1L2TVContext(name, mesh, m; T, tdata, S,
est, g, u, p1, p2, du, dp1, dp2, nablau, nabladu)
end
+function p1_project!(p1, alpha1)
+ p1.space.element == DP0() || p1.space.element == P1() ||
+ p1.space.element == DP1() ||
+ throw(ArgumentError("element unsupported"))
+ p1.data .= clamp.(p1.data, -alpha1, alpha1)
+end
+
+function p2_project!(p2, lambda)
+ p2.space.element::DP1
+ p2d = reshape(p2.data, prod(p2.space.size), :) # no copy
+ for i in axes(p2d, 2)
+ p2in = norm(p2d[:, i])
+ if p2in > lambda
+ p2d[:, i] .*= lambda ./ p2in
+ end
+ end
+end
+
function step!(ctx::L1L2TVContext)
T = ctx.T
S = ctx.S
@@ -151,25 +187,8 @@ function step!(ctx::L1L2TVContext)
ctx.p1.data .+= theta * ctx.dp1.data
ctx.p2.data .+= theta * ctx.dp2.data
- # reproject p1
- function p1_project!(p1, alpha1)
- p1.space.element == DP0() || p1.space.element == P1() ||
- p1.space.element == DP1() ||
- throw(ArgumentError("element unsupported"))
- p1.data .= clamp.(p1.data, -alpha1, alpha1)
- end
+ # reproject p1, p2
p1_project!(ctx.p1, ctx.alpha1)
- # reproject p2
- function p2_project!(p2, lambda)
- p2.space.element::DP1
- p2d = reshape(p2.data, prod(p2.space.size), :) # no copy
- for i in axes(p2d, 2)
- p2in = norm(p2d[:, i])
- if p2in > lambda
- p2d[:, i] .*= lambda ./ p2in
- end
- end
- end
p2_project!(ctx.p2, ctx.lambda)
end
@@ -200,16 +219,6 @@ function step_pd!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
end
interpolate!(p2_next, p2_update; ctx.p2, nablau=nabla(u_new))
- function p2_project!(p2, lambda)
- p2.space.element::DP1
- p2d = reshape(p2.data, prod(p2.space.size), :) # no copy
- for i in axes(p2d, 2)
- p2in = norm(p2d[:, i])
- if p2in > lambda
- p2d[:, i] .*= lambda ./ p2in
- end
- end
- end
p2_project!(p2_next, ctx.lambda)
ctx.dp2.data .= p2_next.data .- ctx.p2.data
@@ -286,25 +295,8 @@ function step_pd2!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
end
interpolate!(p2_next, p2_update; ctx.p2, nablau=nabla(u_new))
- # reproject p1
- function p1_project!(p1, alpha1)
- p1.space.element == DP0() || p1.space.element == P1() ||
- p1.space.element == DP1() ||
- throw(ArgumentError("element unsupported"))
- p1.data .= clamp.(p1.data, -alpha1, alpha1)
- end
+ # reproject p1, p2
p1_project!(p1_next, ctx.alpha1)
- # reproject p2
- function p2_project!(p2, lambda)
- p2.space.element::DP1
- p2d = reshape(p2.data, prod(p2.space.size), :) # no copy
- for i in axes(p2d, 2)
- p2in = norm(p2d[:, i])
- if p2in > lambda
- p2d[:, i] .*= lambda ./ p2in
- end
- end
- end
p2_project!(p2_next, ctx.lambda)
ctx.dp1.data .= p1_next.data .- ctx.p1.data
@@ -450,10 +442,6 @@ function output(ctx::L1L2TVContext, filename, fs...)
return vtk
end
-toimg(arr) = Gray.(clamp.(arr, 0., 1.))
-loadimg(x) = reverse(transpose(Float64.(FileIO.load(x))); dims = 2)
-saveimg(io, x) = FileIO.save(io, transpose(reverse(toimg(x); dims = 2)))
-
function primal_energy(ctx::L1L2TVContext)
function integrand(x; g, u, nablau, tdata)
return ctx.alpha1 * huber(norm(ctx.T(tdata, u) - g), ctx.gamma1) +
@@ -615,14 +603,6 @@ function denoise_pd(img; df=nothing, name, algorithm, params...)
return ctx
end
-"""
- logfilter(dt; a=20)
-
-Reduce dataset such that the resulting dataset would have approximately
-equidistant datapoints on a log-scale. `a` controls density.
-"""
-logfilter(dt; a=20) = filter(:k => k -> round(a*log(k+1)) - round(a*log(k)) > 0, dt)
-
function experiment1(img)
path = "data/pd-comparison"
@@ -688,12 +668,15 @@ function inpaint(img, imgmask; name, params...)
return ctx
end
-function optflow(imgf0, imgf1; name, params...)
+function optflow(ctx)
+ imgf0 = ctx.params.imgf0
+ imgf1 = ctx.params.imgf1
size(imgf0) == size(imgf1) ||
throw(ArgumentError("non-matching dimensions"))
m = 2
- mesh = init_grid(imgf0; type=:vertex)
+ #mesh = init_grid(imgf0; type=:vertex)
+ mesh = init_grid(imgf0, 50, 50)
# optflow specific stuff
Vg = FeSpace(mesh, P1(), (1,))
@@ -705,28 +688,51 @@ function optflow(imgf0, imgf1; name, params...)
T(tdata, u) = tdata * u
S(u, nablau) = nablau
- ctx = L1L2TVContext(name, mesh, m; T, tdata = nablafw, S, params...)
+ name = "run"
+ st = L1L2TVContext(name, mesh, m; T, tdata = nablafw, S,
+ ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta,
+ ctx.params.gamma1, ctx.params.gamma2)
+
+ function reproject!()
+ project_img!(f0, imgf0)
+ project_img!(f1, imgf1)
+ # FIXME: currently dual grid only
+ #interpolate!(f0, x -> imgf0[round.(Int, x)...])
+ #interpolate!(f1, x -> imgf1[round.(Int, x)...])
+ end
+ reproject!()
- # FIXME: currently dual grid only
- interpolate!(f0, x -> imgf0[round.(Int, x)...])
- interpolate!(f1, x -> imgf1[round.(Int, x)...])
fw.data .= f1.data
g_optflow(x; u, f0, fw, nablafw) =
nablafw * u - (fw - f0)
- interpolate!(ctx.g, g_optflow; ctx.u, f0, fw, nablafw)
+ interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw)
save_optflow(i) =
- output(ctx, "output/$(ctx.name)_$(lpad(i, 5, '0')).vtu",
- ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.est, f0, f1, fw)
+ output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
+ st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
- pvd = paraview_collection("output/$(ctx.name).pvd")
+ pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd"))
pvd[0] = save_optflow(0)
for i in 1:10
- step!(ctx)
- estimate!(ctx)
+ step!(st)
+ estimate!(st)
pvd[i] = save_optflow(i)
println()
end
+ vtk_save(pvd)
+
+ #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
+ #saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob))
+ #savedata(ctx, "data.tex"; lambda=λ, beta=β, tau=τ, maxiters, energymin,
+ # width=size(fo, 2), height=size(fo, 1))
+end
+
+function experiment_optflow_middlebury(ctx)
+ imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
+ imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
+
+ ctx = Util.Context(ctx; imgf0, imgf1)
+ optflow(ctx)
return ctx
end
diff --git a/scripts/util.jl b/scripts/util.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f57f0d7b2b58fc1b41f40e173be88292b9dd1d6b
--- /dev/null
+++ b/scripts/util.jl
@@ -0,0 +1,100 @@
+module Util
+
+export savedata
+
+# Path Management
+
+"""
+ Context{P}
+
+Execution context for an experiment.
+"""
+struct Context{P}
+ indir::String
+ outdir::String
+ params::P
+ # TODO: get rid of this when latex output does no longer depend on path
+ path::String
+end
+
+Context(root::String; params...) = Context(root, root, NamedTuple(params), "")
+Context(inroot::String, outroot::String; params...) =
+ Context(inroot, outroot, NamedTuple(params), "")
+Context(ctx::Context, path::String; params...) =
+ Context(joinpath(ctx.indir, path), joinpath(ctx.outdir, path),
+ merge(ctx.params, NamedTuple(params)),
+ joinpath(ctx.path, path))
+Context(ctx::Context; params...) =
+ Context(ctx.indir, ctx.outdir,
+ merge(ctx.params, NamedTuple(params)),
+ ctx.path)
+
+"""
+ (ctx::Context)(f, subdir; params...)
+
+Execute experiment `f` on a context relative to `ctx` and `subdir` using
+parameters `params`.
+"""
+function (ctx::Context)(f::Function, subdir::String; params...)
+ # maybe debug output with params?
+ subctx = Context(ctx, subdir; params...)
+ println("starting experiment `$(subctx.path)` ...")
+ f(subctx)
+end
+
+# LaTeX Data Output
+
+# TODO: don't depend on ctx.path and use filename only
+savedata(ctx::Context, filename; data...) =
+ open(joinpath(ctx.outdir, filename); write=true) do io
+ _savedata(io, ctx.path; data...)
+ end
+
+function _savedata(io, path; data...)
+ isvalidname(x) = contains(string(x), r"^[^/\\]+$")
+
+ # define tex command to access data
+ texcmd = "Data"
+ write(io, """
+ \\providecommand{\\$(texcmd)}[1]{
+ \\csname $(texcmd)/#1\\endcsname
+ }
+ """)
+
+ for (key, value) in pairs(data)
+ isvalidname(key) && isvalidname(value) ||
+ throw(ArgumentError("invalid key/value: $(key => value)"))
+
+ # actual data encoded as tex-commands
+ write(io, "\\expandafter\\def")
+ write(io, "\\csname $(texcmd)/$(path)/$(key)\\endcsname{")
+ show(io, MIME("text/latex"), value)
+ write(io, "}\n")
+ end
+end
+
+# LaTex special printing rules
+
+Base.show(io, ::MIME"text/latex", x) = print(io, x)
+
+function Base.show(io::IO, ::MIME"text/latex", x::AbstractFloat)
+ exponent = x == 0. ? 0 : floor(Int, log10(x))
+ if abs(exponent) <= 3
+ print(io, "\\ensuremath{", x, "}")
+ else
+ base = x / 10. ^ exponent
+ print(io, "\\ensuremath{", base, "\\cdot 10^{", exponent, "}}")
+ end
+end
+
+function Base.show(io::IO, ::MIME"text/latex", x::Integer)
+ p = log10(x)
+ pr = round(Int, p)
+ if x == 10^pr
+ print(io, "\\ensuremath{10^{", pr, "}}")
+ else
+ print(io, "\\ensuremath{", x, "}")
+ end
+end
+
+end
diff --git a/src/SemiSmoothNewton.jl b/src/SemiSmoothNewton.jl
index eae1c6550f9cb6f6be892a1b22f6253291577616..47b41372fbecc66049c15a2184f0bfc212fa6425 100644
--- a/src/SemiSmoothNewton.jl
+++ b/src/SemiSmoothNewton.jl
@@ -5,8 +5,6 @@ include("mesh.jl")
include("function.jl")
include("operator.jl")
-include("run.jl")
-
function clip_line(r0, r1, l0, l1)
d_l = -(l1[1] - l0[1])
d_r = +(l1[1] - l0[1])
diff --git a/src/mesh.jl b/src/mesh.jl
index 86a9f23cfbb72d8b5e72476a91d225983a7616db..12da520bd588f83325710a4a5f5b28387763b08e 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -1,4 +1,7 @@
-export init_grid, init_hgrid, save, refine!, cells, vertices
+export init_grid, init_hgrid, save, refine!, cells, vertices,
+ ndims_domain, ndims_space
+
+using LinearAlgebra: norm
using WriteVTK
using StaticArrays: SVector