The active contour model in the paper “Active Contours Without Edges” by Chan and Vese handles images with statistically different foregrounds and backgrounds by utilizing a region-based approach, instead of relying solely on edge information, which may not be well-defined or present in images with homogeneous regions. This makes it suitable for identifying reflections, lights, or surfaces of water as shown in the demo.

The paper presents a variational formulation for the active contour model, which is based on the level set method. The model is designed to minimize a functional that consists of two terms: a data term that measures the fit between the image and the evolving curve, and a regularization term that encourages smoothness of the curve. The data term is based on the intensity values of the image, while the regularization term is based on the length of the curve and the area enclosed by it. The paper also presents an efficient numerical scheme for solving the variational formulation, which is based on the level set method and the gradient descent method. The proposed model is tested on a variety of images, including synthetic and real images, and is shown to be effective in detecting and segmenting objects in images with complex boundaries.

The algorithm works by evolving a curve that separates the foreground and background regions of the image. The curve is represented as a signed distance function (SDF) `phi`

, which is initialized from the binary mask `init_mask`

. At each iteration, the algorithm updates the curve by minimizing an energy functional that consists of two terms: a data term that measures the fit between the image and the curve, and a regularization term that encourages smoothness of the curve. The data term is based on the difference between the mean intensities of the foreground and background regions, while the regularization term is based on the curvature of the curve. The energy functional is minimized using a gradient descent method, and the curve is evolved by updating the SDF using the computed gradient. The algorithm also maintains the Courant-Friedrichs-Lewy (CFL) condition to ensure stability of the numerical scheme.

So what does all this look like? Well the code ends up being fairly verbose, so let me summarize the overall algorithm.

```
function seg = region_seg(I,init_mask,max_its,alpha,display)
if alpha is not given, set alpha to 0.2
if display is not given, set display to true
Convert I to a 2D double matrix
Create a signed distance map (SDF) from init_mask using mask2phi
for its = 1 to max_its
Get the curve's narrow band by finding phi values between -1.2 and 1.2
Find the interior and exterior mean values u and v
Calculate the force F from image information
Calculate the curvature penalty
Calculate the gradient descent dphidt to minimize energy
Calculate the time step dt to maintain the CFL condition
Evolve the curve by updating phi using dphidt and dt
Keep the SDF smooth using sussman
If display is true and its is a multiple of 20, show the intermediate output
end
If display is true, show the final output
Create a mask from the SDF by setting seg to phi<=0
Return seg
end
```

And here is a C++ Unit Test using OpenCV

```
#include <iostream>
#include <opencv2/opencv.hpp>
cv::Mat region_seg(cv::Mat I, cv::Mat init_mask, int max_its, double alpha = 0.2, bool display = true) {
// Convert I to a 2D double matrix
cv::Mat I_double;
I.convertTo(I_double, CV_64F);
// Create a signed distance map (SDF) from init_mask using mask2phi
cv::Mat phi = mask2phi(init_mask);
for (int its = 1; its <= max_its; its++) {
// Get the curve's narrow band by finding phi values between -1.2 and 1.2
cv::Mat narrow_band = get_narrow_band(phi, -1.2, 1.2);
// Find the interior and exterior mean values u and v
double u = calculate_mean(narrow_band, I_double, true);
double v = calculate_mean(narrow_band, I_double, false);
// Calculate the force F from image information
cv::Mat F = calculate_force(I_double, u, v, alpha);
// Calculate the curvature penalty
cv::Mat curvature_penalty = calculate_curvature_penalty(phi);
// Calculate the gradient descent dphidt to minimize energy
cv::Mat dphidt = calculate_gradient_descent(F, curvature_penalty);
// Calculate the time step dt to maintain the CFL condition
double dt = calculate_time_step(phi, dphidt);
// Evolve the curve by updating phi using dphidt and dt
phi += dt * dphidt;
// Keep the SDF smooth using sussman
sussman(phi);
// If display is true and its is a multiple of 20, show the intermediate output
if (display && its % 20 == 0) {
cv::imshow("Intermediate Output", phi);
cv::waitKey(0);
}
}
// If display is true, show the final output
if (display) {
cv::imshow("Final Output", phi);
cv::waitKey(0);
}
// Create a mask from the SDF by setting seg to phi<=0
cv::Mat seg = phi <= 0;
return seg;
}
int main() {
// Example usage
cv::Mat I = cv::imread("image.jpg", cv::IMREAD_GRAYSCALE);
cv::Mat init_mask = cv::imread("mask.jpg", cv::IMREAD_GRAYSCALE);
cv::Mat seg = region_seg(I, init_mask, 100);
return 0;
}
```