|
29 | 29 | "- $f = m u_{\\mathrm{true}}$: observed incomplete image\n", |
30 | 30 | "- $u$: reconstructed image\n", |
31 | 31 | "\n", |
32 | | - "We compute $u$ by minimising\n", |
| 32 | + "We compute $u$ by minimizing\n", |
33 | 33 | "\n", |
34 | 34 | "$$\n", |
35 | 35 | "J(u)= {1 \\over 2}\\beta \\int_\\Omega m(u-f)^2\\,\\mathrm{d}x\n", |
36 | 36 | "+ \\alpha \\int_\\Omega \\sqrt{||\\nabla u||^2 + \\varepsilon^2}~\\mathrm{d}x.\n", |
37 | 37 | "$$\n", |
38 | 38 | "\n", |
39 | 39 | "The first term enforces agreement with the known image data, while\n", |
40 | | - "the second term is a smoothed total variation regularisation term.\n", |
| 40 | + "the second term is a smoothed total variation regularization term.\n", |
41 | 41 | "It promotes piecewise smooth solution and preserves edges\n", |
42 | 42 | "$\\alpha$ and $\\beta$ control the balance between the data fidelity\n", |
43 | 43 | "(fit to f) and smoothness.\n", |
|
106 | 106 | "source": [ |
107 | 107 | "We use first order Lagrange elements for discretizing the image.\n", |
108 | 108 | "In this space, the DOFs are the values of u at mesh vertices\n", |
109 | | - "the solution is continous but has piecewise constant gradient" |
| 109 | + "the solution is continuous but has piecewise constant gradient" |
110 | 110 | ] |
111 | 111 | }, |
112 | 112 | { |
|
175 | 175 | "\n", |
176 | 176 | "We construct a mask with random \"holes\" inside the square\n", |
177 | 177 | "* small circular regions are removed and set to 0\n", |
178 | | - "* everyhere else remains known (1)" |
| 178 | + "* everywhere else remains known (1)" |
179 | 179 | ] |
180 | 180 | }, |
181 | 181 | { |
|
209 | 209 | " # number of speckles\n", |
210 | 210 | " num_speckles = 25\n", |
211 | 211 | " # random centers\n", |
212 | | - " generator = np.random.Generator(np.random.MT19937(0)) # random seed for reproducibility\n", |
| 212 | + " generator = np.random.Generator(\n", |
| 213 | + " np.random.MT19937(0)\n", |
| 214 | + " ) # random seed for reproducibility\n", |
213 | 215 | "\n", |
214 | 216 | " cx = generator.uniform(0.25, 0.75, num_speckles)\n", |
215 | 217 | " cy = generator.uniform(0.25, 0.75, num_speckles)\n", |
216 | 218 | " # random radii (small + varied)\n", |
217 | | - " radii = generator.uniform(0.012, 0.035, num_speckles)\n", |
| 219 | + " radii = generator.uniform(0.012, 0.035, num_speckles)\n", |
218 | 220 | " # create holes. mask =0 inside circles\n", |
219 | 221 | " for i in range(num_speckles):\n", |
220 | 222 | " r2 = (X - cx[i]) ** 2 + (Y - cy[i]) ** 2\n", |
|
281 | 283 | "\n", |
282 | 284 | "where # $\\varepsilon$ is the smoothing of the TV:\n", |
283 | 285 | "* large $\\varepsilon$ smoother more like quadratic diffusion\n", |
284 | | - "* small $\\varepsilon$ closer to true TV edge pereserving" |
| 286 | + "* small $\\varepsilon$ closer to true TV edge preserving" |
285 | 287 | ] |
286 | 288 | }, |
287 | 289 | { |
|
302 | 304 | "metadata": {}, |
303 | 305 | "source": [ |
304 | 306 | "Smoothed TV inpainting energy functional.\n", |
305 | | - "We define the energy J(u) and use ufl.derivative to obtain\n", |
306 | | - "the residual form F.\n", |
| 307 | + "We define the energy J(u) and use {py:func}`ufl.derivative` to obtain\n", |
| 308 | + "the residual form $F(u; v)=0 \\quad\\forall v\\in V$:\n", |
307 | 309 | "\n", |
308 | 310 | "$$\n", |
309 | 311 | "J(u) = {1 \\over 2}\\beta\\int_\\Omega m(u-f)^2\\,dx\n", |
310 | | - "+ \\alpha\\int_\\Omega \\sqrt{||\\nabla u||^2+\\varepsilon^2}\\,dx\n", |
| 312 | + "+ \\alpha\\int_\\Omega \\sqrt{||\\nabla u||^2+\\varepsilon^2}~\\mathrm{d}x\n", |
311 | 313 | "$$\n", |
312 | 314 | "\n", |
313 | 315 | "Taking the first variation gives the weak form F(u; v).\n", |
|
318 | 320 | "+ \\alpha\\int_\\Omega\n", |
319 | 321 | "{\\nabla u\\cdot\\nabla v\n", |
320 | 322 | "\\over\n", |
321 | | - "\\sqrt{||\\nabla u||^2+\\varepsilon^2}}\\,dx\n", |
| 323 | + "\\sqrt{||\\nabla u||^2+\\varepsilon^2}}~\\mathrm{d}x\n", |
322 | 324 | "= 0 \\quad \\forall v\\in V.\n", |
323 | 325 | "$$" |
324 | 326 | ] |
|
330 | 332 | "metadata": {}, |
331 | 333 | "outputs": [], |
332 | 334 | "source": [ |
333 | | - "v = ufl.TestFunction(V)" |
334 | | - ] |
335 | | - }, |
336 | | - { |
337 | | - "cell_type": "code", |
338 | | - "execution_count": null, |
339 | | - "id": "0055665b", |
340 | | - "metadata": {}, |
341 | | - "outputs": [], |
342 | | - "source": [ |
| 335 | + "v = ufl.TestFunction(V)\n", |
| 336 | + "\n", |
343 | 337 | "J_energy = (\n", |
344 | 338 | " 0.5 * beta * m * (u - f) ** 2 * ufl.dx\n", |
345 | 339 | " + alpha * ufl.sqrt(ufl.inner(ufl.grad(u), ufl.grad(u)) + eps**2) * ufl.dx\n", |
346 | | - ")" |
347 | | - ] |
348 | | - }, |
349 | | - { |
350 | | - "cell_type": "code", |
351 | | - "execution_count": null, |
352 | | - "id": "3fbade6e", |
353 | | - "metadata": {}, |
354 | | - "outputs": [], |
355 | | - "source": [ |
| 340 | + ")\n", |
| 341 | + "\n", |
356 | 342 | "F = ufl.derivative(J_energy, u, v)" |
357 | 343 | ] |
358 | 344 | }, |
359 | 345 | { |
360 | 346 | "cell_type": "markdown", |
361 | | - "id": "77b7387d", |
| 347 | + "id": "163eb995", |
362 | 348 | "metadata": {}, |
363 | 349 | "source": [ |
364 | | - "This formulation is based on total variation (TV) regulaization\n", |
| 350 | + "This formulation is based on total variation (TV) regularization\n", |
365 | 351 | "for image denoising and inpainting\n", |
366 | 352 | "{cite:t}`tv-RUDIN1992TV,tv-CHAN2001TV`." |
367 | 353 | ] |
|
423 | 409 | "source": [ |
424 | 410 | "FEM Metrics\n", |
425 | 411 | "Global number of degrees of freedom reports the size of the\n", |
426 | | - "finite element discretisation H1 seminorm error measures the\n", |
| 412 | + "finite element discretization H1 seminorm error measures the\n", |
427 | 413 | "gradient error\n", |
428 | 414 | "\n", |
429 | 415 | "$$\n", |
430 | 416 | " \\vert\\vert\\nabla(u-u_{true})\\vert\\vert_{L_2 (\\Omega)}\n", |
431 | 417 | "$$\n", |
432 | 418 | "\n", |
433 | | - "This is useful as TV regulization is gradient based.\n", |
| 419 | + "This is useful as TV regularization is gradient based.\n", |
434 | 420 | "Smaller values mean the reconstruction recovers edge structure better" |
435 | 421 | ] |
436 | 422 | }, |
|
622 | 608 | "$$\n", |
623 | 609 | "\n", |
624 | 610 | "A decrease in the objective show that the nonlinear optimization\n", |
625 | | - "improved the damaged image undeer the smoothed TV model" |
| 611 | + "improved the damaged image under the smoothed TV model" |
626 | 612 | ] |
627 | 613 | }, |
628 | 614 | { |
|
654 | 640 | "J0 = 0.5 * float(beta) * J0_data + float(alpha) * J0_tv" |
655 | 641 | ] |
656 | 642 | }, |
| 643 | + { |
| 644 | + "cell_type": "markdown", |
| 645 | + "id": "f01dfcb3", |
| 646 | + "metadata": {}, |
| 647 | + "source": [ |
| 648 | + "Printing statements for validation and metrics\n", |
| 649 | + "If on main process" |
| 650 | + ] |
| 651 | + }, |
657 | 652 | { |
658 | 653 | "cell_type": "code", |
659 | 654 | "execution_count": null, |
660 | 655 | "id": "65e16ee1", |
661 | 656 | "metadata": {}, |
662 | 657 | "outputs": [], |
663 | 658 | "source": [ |
664 | | - "# Printing statments for validation and metrics\n", |
665 | | - "# If on main process\n", |
666 | 659 | "if msh.comm.rank == 0:\n", |
667 | 660 | " print(\"---Smoothed TV inpainting results---\")\n", |
668 | 661 | "\n", |
|
699 | 692 | "We construct fields that allow us to visually asses the quality\n", |
700 | 693 | "of the reconstruction\n", |
701 | 694 | "$u-u_{true}$ is the global reconstruction error\n", |
702 | | - "$(1-m)(u-u_{true})$ is the hole error, restriced to the missing regions" |
| 695 | + "$(1-m)(u-u_{true})$ is the hole error, restricted to the missing regions" |
703 | 696 | ] |
704 | 697 | }, |
705 | 698 | { |
|
727 | 720 | "To plot in matplotlib\n", |
728 | 721 | "1. extract the coordinates of the DOFs\n", |
729 | 722 | "2. extract the function-space dofmap connectivity (triangles)\n", |
730 | | - "3. build a Triangulation object\n", |
| 723 | + "3. build a {py:class}`Triangulation<matplotlib.tri.Triangulation>` object\n", |
| 724 | + "\n", |
731 | 725 | "This allows matplotlib to render the piecewise linear FEM solution" |
732 | 726 | ] |
733 | 727 | }, |
|
739 | 733 | "outputs": [], |
740 | 734 | "source": [ |
741 | 735 | "coords = V.tabulate_dof_coordinates()\n", |
742 | | - "x, y = coords[:, 0], coords[:, 1]" |
743 | | - ] |
744 | | - }, |
745 | | - { |
746 | | - "cell_type": "code", |
747 | | - "execution_count": null, |
748 | | - "id": "8a5a2ec4", |
749 | | - "metadata": { |
750 | | - "lines_to_next_cell": 2 |
751 | | - }, |
752 | | - "outputs": [], |
753 | | - "source": [ |
| 736 | + "x, y = coords[:, 0], coords[:, 1]\n", |
| 737 | + "\n", |
754 | 738 | "triangles = V.dofmap.list\n", |
755 | 739 | "triang = mtri.Triangulation(x, y, triangles)" |
756 | 740 | ] |
|
761 | 745 | "metadata": {}, |
762 | 746 | "source": [ |
763 | 747 | "Plotting\n", |
764 | | - "We use tripcolor to plot scalar fields defined on a triangulated mesh\n", |
765 | | - "shading= \"flat\" shows piecewise constant coloring per triangle\n", |
| 748 | + "We use {py:func}`matplotlib.pyplot.tripcolor` to plot scalar fields defined on a\n", |
| 749 | + "triangulated mesh shading= \"flat\" shows piecewise constant coloring per triangle\n", |
766 | 750 | "which better reflects the discrete FEM representations" |
767 | 751 | ] |
768 | 752 | }, |
|
796 | 780 | "lim = np.max(np.abs(u_minus_u_true.x.array))\n", |
797 | 781 | "# $u-u_{true}$ is the global reconstruction error\n", |
798 | 782 | "plot_field(\n", |
799 | | - " axes[1, 1], u_minus_u_true.x.array, \"u - u_true\", fig, cmap=\"coolwarm\", vmin=-lim, vmax=lim\n", |
| 783 | + " axes[1, 1],\n", |
| 784 | + " u_minus_u_true.x.array,\n", |
| 785 | + " \"u - u_true\",\n", |
| 786 | + " fig,\n", |
| 787 | + " cmap=\"coolwarm\",\n", |
| 788 | + " vmin=-lim,\n", |
| 789 | + " vmax=lim,\n", |
800 | 790 | ")\n", |
801 | 791 | "# Hole only errors\n", |
802 | 792 | "lim = np.max(np.abs(hole_error_field.x.array))\n", |
|
830 | 820 | ], |
831 | 821 | "metadata": { |
832 | 822 | "jupytext": { |
| 823 | + "formats": "ipynb,py:light", |
833 | 824 | "main_language": "python" |
834 | 825 | } |
835 | 826 | }, |
|
0 commit comments