Minimization in Chan-Vese convex relation model

I am working on an assignment with the Chan-Vese convex relaxation model for two-phase image segmentation, i.e., coming up with two regions — the object and the background. Specifically, I need to implement an algorithm that minimizes the energy function

$$\min_{\theta \in [0,1], c_1, c_2} \int|\nabla\theta| + \int_\Omega(u_0(x,y)-c_1)^2\theta(x,y)dxdy + \int_\Omega(u_0(x,y)-c_2)^2(1-\theta(x,y))dxdy$$

where

• $$u_0$$ is the grayscale image

• $$c_1$$ is the average intensity inside the boundary (the object)

• $$c_2$$ is the average intensity outside (the background)

• $$\Omega$$ is the domain

• $$\theta(x,y)$$ is a value between 0 and 1 which ultimately will give the object and the background. If $$\theta \geq \frac 12$$ we will call that the object and $$\theta \leq \frac 12$$ will be the background.

My problem lies in figuring out an updating scheme for $$\theta$$. The general approach I've learned so far is to differentiate with respect to which every variable, set the expression equal to zero, and that gives the updating rule. For example, differentiating wrt $$c_1$$ and setting equal to zero gives

$$c_1 = \frac{\int_\Omega u_0 \theta}{\int_\Omega \theta}$$

However, I am lost on how to update $$\theta$$. We have done lots of minimization of expressions like $$\|Ax-b\|_2^2$$ but here I don't know how to turn the integral

$$\int_\Omega(u_0(x,y)-c_1)^2\theta(x,y)dxdy$$

into some norm squared. Without the $$\theta$$ I know the integral is the $$L^2$$ norm for functions, i.e.,

$$\int_\Omega(u_0(x,y)-c_1)^2dxdy = \|u_0-c_1\|_2^2$$

As a suggestion I was told to think of $$\theta$$ operating pointwise but I don't know how that helps.