Manipulating images

Hi everyone,
I am trying to solve the heat equation using a matrix u0 of pixel values as initial condition. I am struggling with the manipulation of images.
Is there any complete template code to read, store, …show an image and define it as initial condition in Fenics?
Many thinks for your help.

1 Like

Hi, I used something like:

import matplotlib.pyplot as plt
img = plt.imread("my_image.png")
(Nx, Ny) = img.shape

mesh = UnitSquareMesh(Nx, Ny, "crossed")

class FE_image(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        p = Cell(mesh, ufc_cell.index).midpoint()
        i, j = int(p[0]*(Nx-1)), int(p[1]*(Ny-1))
        value[:] = img[-(j+1), i]

    def value_shape(self):
        return ()

y = FE_image()

y is then an Expression that you can interpolate or project as you want.

1 Like

Thank you very much for your reply. I tried to apply your code but still don’t understand how to use the ‘y’ expression in Fenics. Any more sights please…

Without any more precision of what you really want to do we cannot really help you precisely. As mentionned you can interpolate it on a FunctionSpace V to use it as initial condition:

u = Function(V)
u.interpolate(y)

I tried the following but I did not understand why I got the value 30 for all the vertices.
My goal is to solve the time-dependent heat equation using an image as initial condition.

from fenics import *
import matplotlib.pyplot as plt
import numpy as np
img = np.array([[0, 10], [30, 90]])
(Nx, Ny) = img.shape
mesh = UnitSquareMesh(Nx, Ny, “crossed”)

class FE_image(UserExpression):
def eval_cell(self, value, x, ufc_cell):
p = Cell(mesh, ufc_cell.index).midpoint()
i, j = int(p[0](Nx-1)), int(p[1](Ny-1))
value[:] = img[-(j+1), i]

def value_shape(self):
    return ()

y = FE_image()

V = FunctionSpace(mesh, ‘P’, 1)
#V = FunctionSpace(mesh, “Lagrange”, 1)

u0 = Function(V)
u0.interpolate(y)
#u0 = project(y, V)
#u0 = interpolate(y, V)
#u0 = interpolate(Expression(‘x[0] + x[1]’, degree=1), V)

vertex_values = u0.compute_vertex_values()
coordinates = mesh.coordinates()

for i, x in enumerate(coordinates):
… print(‘vertex %d: vertex_values[%d] = %g\tu(%s) = %g’ %
… (i, i, vertex_values[i], x, u0(x)))

vertex 0: vertex_values[0] = 30 u([0. 0.]) = 30
vertex 1: vertex_values[1] = 30 u([0.5 0. ]) = 30
vertex 2: vertex_values[2] = 30 u([1. 0.]) = 30
vertex 3: vertex_values[3] = 30 u([0. 0.5]) = 30
vertex 4: vertex_values[4] = 30 u([0.5 0.5]) = 30
vertex 5: vertex_values[5] = 30 u([1. 0.5]) = 30
vertex 6: vertex_values[6] = 30 u([0. 1.]) = 30
vertex 7: vertex_values[7] = 30 u([0.5 1. ]) = 30
vertex 8: vertex_values[8] = 30 u([1. 1.]) = 30
vertex 9: vertex_values[9] = 30 u([0.25 0.25]) = 30
vertex 10: vertex_values[10] = 30 u([0.75 0.25]) = 30
vertex 11: vertex_values[11] = 30 u([0.25 0.75]) = 30
vertex 12: vertex_values[12] = 30 u([0.75 0.75]) = 30

1 Like

Sorry there was a mistake in the FE_image class, change to:

i, j = int(p[0]*Nx), int(p[1]*Ny)
1 Like

Much better, thank you.
Is it possible to change the code to get structured data?
For example, if I have a 3 by 3 image, I want to get 9 vertex_values.
I want to access a value by its i and j indices, i counting cells in the x direction, and j counting cells in the y direction.

Thanks for this, it is very helpful. Just a note: is it working only for square images? If Nx=Ny it works fine, otherwise I get an out-of-bounds error. How does this sound to you?

Thanks

Yes sorry, you might have to sxchange Nx and Ny in the definition of i and j

1 Like