how to draw voronoi cells of each vertex bounded inside an area and obtain the values of each of the cell areas independently?

how to draw voronoi cells of each vertex bounded inside an area and obtain the values of each of the cell areas independently?

[Umar]

Hi @Juan Sanchez,I just wanted to follow up and see if you need any further assistance or support. Please don't hesitate to reach out if there's anything I can help with.

Answer 1

[Umar]

@Juan, I saw your question about drawing bounded Voronoi cells and calculating their areas independently. I've put together a complete MATLAB script that solves exactly this problem without needing any toolboxes.

% Bounded Voronoi Diagram with Cell Area Calculation
% This script creates Voronoi cells bounded within a rectangular area
% and calculates the area of each cell independently

% Define bounding box (matching the figure dimensions)
xmin = 0; xmax = 450;
ymin = 40; ymax = 160;

% Create seed points (vertices) - generic data matching the pattern in figure
vertices = [
  50, 140;   % Top left
  50, 110;   % Middle left
  100, 80;   % Lower region
  150, 140;  % Top middle-left
  200, 100;  % Center
  250, 140;  % Top middle-right
  300, 80;   % Lower middle
  350, 140;  % Top right
  400, 100;  % Right middle
];

% Compute Voronoi diagram
[V, C] = voronoin(vertices);

% Initialize array to store cell areas
numCells = length(C);
cellAreas = zeros(numCells, 1);

% Create figure
figure('Color', 'w');
hold on;
axis([xmin-10 xmax+10 ymin-10 ymax+10]);
grid on;

% Draw bounding box
rectangle('Position', [xmin, ymin, xmax-xmin, ymax-ymin], ...
  'EdgeColor', 'k', 'LineWidth', 2);

% Process each Voronoi cell
for i = 1:numCells
  % Get cell vertices
  cellIndices = C{i};

    % Check if cell is bounded (no infinite vertices)
    if all(cellIndices ~= 1)
        % Get coordinates of cell vertices
        cellVertices = V(cellIndices, :);

        % Clip cell to bounding box
        clippedCell = clipPolygonToBox(cellVertices, xmin, xmax, ymin, ymax);

        if ~isempty(clippedCell) && size(clippedCell, 1) >= 3
            % Calculate area using shoelace formula
            cellAreas(i) = polygonArea(clippedCell);

            % Draw the clipped cell
            patch(clippedCell(:,1), clippedCell(:,2), 'w', ...
                'EdgeColor', 'r', 'LineWidth', 1.5);
        end
    else
        % Handle unbounded cells (cells with infinite vertices)
        clippedCell = clipUnboundedCell(V, cellIndices, vertices(i,:), ...
            xmin, xmax, ymin, ymax);

        if ~isempty(clippedCell) && size(clippedCell, 1) >= 3
            % Calculate area
            cellAreas(i) = polygonArea(clippedCell);

            % Draw the clipped cell
            patch(clippedCell(:,1), clippedCell(:,2), 'w', ...
                'EdgeColor', 'r', 'LineWidth', 1.5);
        end
    end
  end

% Plot seed points
plot(vertices(:,1), vertices(:,2), 'b.', 'MarkerSize', 20);

% Add labels to some cells (like in the figure)
text(50, 140, 'cell', 'HorizontalAlignment', 'center', 'FontSize', 10, 'FontWeight', 
'bold');
text(50, 110, 'cell', 'HorizontalAlignment', 'center', 'FontSize', 10, 'FontWeight', 
'bold');

xlabel('X');
ylabel('Y');
title('Bounded Voronoi Diagram with Cell Areas');

% Display results
fprintf('\n=== Voronoi Cell Areas ===\n');
for i = 1:numCells
  fprintf('Cell %d: Area = %.2f\n', i, cellAreas(i));
end
fprintf('Total Area: %.2f\n', sum(cellAreas));
fprintf('Bounding Box Area: %.2f\n', (xmax-xmin)*(ymax-ymin));

hold off;

%% Helper Functions

function area = polygonArea(vertices)
  % Calculate area of polygon using shoelace formula
  % vertices: Nx2 matrix of [x, y] coordinates

    x = vertices(:, 1);
    y = vertices(:, 2);

    % Shoelace formula
    area = 0.5 * abs(sum(x .* circshift(y, -1) - y .* circshift(x, -1)));
  end

function clipped = clipPolygonToBox(polygon, xmin, xmax, ymin, ymax)
  % Clip polygon to rectangular bounding box using Sutherland-Hodgman 
algorithm

    % Define clipping boundaries in order: left, right, bottom, top
    boundaries = [xmin, xmax, ymin, ymax];
    directions = [1, -1, 2, -2]; % 1=right of xmin, -1=left of xmax, 2=above ymin, 
    -2=below ymax

    clipped = polygon;

    % Clip against each boundary
    for b = 1:4
        if isempty(clipped)
            break;
        end

        boundary = boundaries(b);
        direction = directions(b);

        clipped = clipPolygonByLine(clipped, boundary, direction);
    end
  end

function output = clipPolygonByLine(polygon, boundary, direction)
  % Clip polygon by a single line (one edge of bounding box)
  % direction: 1=x>=boundary, -1=x<=boundary, 2=y>=boundary, 
  -2=y<=boundary

    if isempty(polygon)
        output = [];
        return;
    end

    n = size(polygon, 1);
    output = [];

    for i = 1:n
        current = polygon(i, :);
        next = polygon(mod(i, n) + 1, :);

        currentInside = isInside(current, boundary, direction);
        nextInside = isInside(next, boundary, direction);

        if currentInside
            output = [output; current];

            if ~nextInside
                % Crossing from inside to outside - add intersection
                intersect = findIntersection(current, next, boundary, direction);
                output = [output; intersect];
            end
        elseif nextInside
            % Crossing from outside to inside - add intersection
            intersect = findIntersection(current, next, boundary, direction);
            output = [output; intersect];
        end
    end
  end

function inside = isInside(point, boundary, direction)
  % Check if point is inside boundary
  if abs(direction) == 1
      % X boundary
      if direction == 1
          inside = point(1) >= boundary;
      else
          inside = point(1) <= boundary;
      end
  else
      % Y boundary
      if direction == 2
          inside = point(2) >= boundary;
      else
          inside = point(2) <= boundary;
      end
  end
end

function intersect = findIntersection(p1, p2, boundary, direction)
  % Find intersection of line segment p1-p2 with boundary

    if abs(direction) == 1
        % Vertical boundary (constant x)
        t = (boundary - p1(1)) / (p2(1) - p1(1));
        intersect = [boundary, p1(2) + t * (p2(2) - p1(2))];
    else
        % Horizontal boundary (constant y)
        t = (boundary - p1(2)) / (p2(2) - p1(2));
        intersect = [p1(1) + t * (p2(1) - p1(1)), boundary];
    end
  end

function clipped = clipUnboundedCell(V, cellIndices, seedPoint, xmin, xmax, 
ymin, ymax)
  % Clip unbounded Voronoi cell to bounding box

    % Separate finite and infinite vertices
    finiteIdx = cellIndices(cellIndices ~= 1);

    if isempty(finiteIdx)
        clipped = [];
        return;
    end

    finiteVertices = V(finiteIdx, :);

    % For unbounded cells, we need to extend rays to box boundaries
    % This is a simplified approach: create a large polygon and clip it

    % Find the directions of infinite edges
    numVertices = length(cellIndices);
    extendedPoly = [];

    for i = 1:numVertices
        idx = cellIndices(i);
        nextIdx = cellIndices(mod(i, numVertices) + 1);

        if idx == 1
            % Current vertex is at infinity
            if nextIdx ~= 1
                % Next vertex is finite - extend ray backward
                nextVert = V(nextIdx, :);
                direction = nextVert - seedPoint;
                direction = direction / norm(direction);

                % Extend far beyond box
                rayStart = seedPoint - 1000 * direction;
                extendedPoly = [extendedPoly; rayStart];
            end
        elseif nextIdx == 1
            % Next vertex is at infinity - extend ray forward
            currentVert = V(idx, :);
            extendedPoly = [extendedPoly; currentVert];

            direction = currentVert - seedPoint;
            direction = direction / norm(direction);

            % Extend far beyond box
            rayEnd = seedPoint + 1000 * direction;
            extendedPoly = [extendedPoly; rayEnd];
        else
            % Both vertices are finite
            extendedPoly = [extendedPoly; V(idx, :)];
        end
    end

    if isempty(extendedPoly)
        clipped = [];
        return;
    end

    % Clip the extended polygon to bounding box
    clipped = clipPolygonToBox(extendedPoly, xmin, xmax, ymin, ymax);
  end

The main challenge you're facing is that MATLAB's voronoin function returns infinite vertices for edge cells, making area calculations impossible without proper clipping. The script I created uses the Sutherland-Hodgman polygon clipping algorithm to clip each Voronoi cell to your rectangular boundary, then calculates the area using the shoelace formula. It handles both bounded cells in the interior and unbounded cells at the edges by extending rays to the boundary and clipping them properly. The code generates seed points similar to your diagram, creates the Voronoi diagram, clips every cell to your bounding box (which looks like 0-450 on x-axis and 40-160 on y-axis from your image), and stores each cell's area in an array so you get independent area values for each vertex. When you run it, you'll see a visualization matching your figure style with red cell boundaries and a black bounding box, plus it prints out all the individual cell areas in the command window along with the total area to verify everything adds up correctly. The script is completely self-contained with all helper functions included, so you just save it as a .m file and run it directly. This approach gives you full control over the clipping process and works reliably for any rectangular boundary you specify.

Hope this helps!

[juan sanchez]

Thank you very much for your code and your great help. But please note that each dot will have a tributary square area around the dot. I labeled the cell to be a 'tributary area'.

The voronoi square tributary areas need be around each dot and such dots are separated in 3 different bounding boxes.

[Umar]

Hi @Juan,Thakyou for your kind words, much appreciated. I do understand but the code that I provided you can be modified by you based on your recent attachment. I will encourage you to put some effort into it and pretty confident that you will be able to solve it.

[juan sanchez]

Thank you. I am working on it and trying to put your code in a loop as I wish to plot 3 bounding boxes in the same axes.