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;
}