how to use fenics image processing because image is an 2D matrix i.e we don't have an expression for f (Poisson eqn) ?

Answer 1

Hi, consider adapting this example

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

N = 40
mesh = UnitSquareMesh(N, N)
x = mesh.coordinates().reshape((-1, 2))

# Vertex with cordinates x[i] has a value in the image at ii[i], jj[i]
h = 1./N
ii, jj = x[:, 0]/h, x[:, 1]/h
ii = np.array(ii, dtype=int)
jj = np.array(jj, dtype=int)

# Image - n by n grid
n = int(sqrt(mesh.num_vertices()))
X, Y = np.meshgrid(np.linspace(0, 1, n), np.linspace(0, 1, n))
image = np.abs(X-Y-1)

# Turn image into CG1 function
# Values are vertex ordered here
image_values = image[ii, jj] 
V = FunctionSpace(mesh, 'CG', 1)
image_f = Function(V)
# Values will be dof ordered
d2v = dof_to_vertex_map(V)
image_values = image_values[d2v]
image_f.vector()[:] = image_values

# Image manip
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = inner(image_f, v)*dx
# Setup solver
A = assemble(a)
solver = KrylovSolver(A, 'cg', 'amg')
# I don't know bcs so the kernel removed at LA level
z = Vector(image_f.vector())
V.dofmap().set(z, 1.0)
z *= z/z.norm('l2')
Z = VectorSpaceBasis([z])
solver.set_nullspace(Z)

# Smooth
for i in range(4):
    b = assemble(L)
    Z.orthogonalize(b)
    solver.solve(image_f.vector(), b)

# Get back to 'image'
v2d = vertex_to_dof_map(V)
new_image = np.zeros_like(image)
values = image_f.vector().array()[v2d]

for (i, j, v) in zip(ii, jj, values): new_image[i, j] = v

plt.figure()
plt.contour(image)

plt.figure()
plt.contour(new_image)

plt.show() 

Comments

ImportError: No module named matplotlib.pyplot, This error is coming then wht can I do?

---

Remove the import statement for matplotlib and the part of code after plt.figure(). The plotting code was there to show that after smoothing the contours went from straight lines to some curves.

---

what is ii and jj?

---

If I give 'DG' instead of 'CG' will it work or somw modification is needed?

---

It depends on the polynomial degree. For piecewise constant fields you can manually build mappings which provide analogous functionality to vertex_to_dof_map and dof_to_vertex_map (think midpoints). For higher degrees, same vertex is shared by a number of dofs, so the logic of this map is not straightforward.

---

Can you help me with the code for Dg for polynomail degree 2 for the same problem.