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.

[juan sanchez]

Hi Umar, Thank you for checking up on this question. I was able to create the voronoi triangles inside a loop for an example of 3 bounding boxes. I used your code but with the dots re-arranged in a rectangular fashion. However, I am only getting irregular polygons for voronoi diagrams instead of rectangular voronoi diagrams. The Voronoi diagrams I am looking to do should be rectangular around each dot but I can not get any of these functions to plot the voronoi diagrams as rectangular. Do you know how can we enforce the voronoi diagrams to be rectangular when the dots are rectangularly positioned? Any information will greatly help. Thank you.

[Umar]

Hi @Juan,

I reviewed your question about getting rectangular Voronoi diagrams when your seed points are arranged in a rectangular grid. I need to clarify something important: standard Voronoi diagrams, including those created by MATLAB's voronoin function that my code uses, mathematically cannot produce rectangular cells even when points are arranged rectangularly. Voronoi cells are defined by perpendicular bisectors between neighboring points, which naturally creates polygonal regions like hexagons or irregular polygons, not rectangles. The code I provided is working correctly and producing the mathematically correct Voronoi diagram for your points. If you truly need rectangular cells around each point, you're not looking for a Voronoi diagram at all but rather a simple grid tessellation where you just divide your bounding box into equal rectangular regions. Alternatively, if you want something closer to rectangular shapes, you could explore Manhattan distance Voronoi diagrams which use city-block distance instead of straight-line distance and produce more square-like growing regions, though MATLAB doesn't have a built-in function for this and you'd need to implement it from scratch. Could you clarify what your actual goal is so I can point you toward the right solution, whether that's a grid tessellation, a different distance metric, or something else entirely?

[juan sanchez]

Thank you for your answer. Looking for information I was able to find this code VoronoiLimitRectSquare at https://nl.mathworks.com/matlabcentral/fileexchange/106380-voronoi-bounded-rectangular-or-square-limits which allows the voronoi to be drawn as rectangles. However this code only extends the lines similarly to your code to the edges of the bounding box. I was able to do the rectangular plot but the problem is that I can not seem to find a method to calculate the areas of each rectangular cell. Which method would you use to find the area using this code? and what difference do you see between VoronoiLimitRectSquare and your provided code? Thank you for your help.

[Umar]

Hi @Juan,

Thanks for your message. One thing to keep in mind is that the VoronoiLimitRectSquare function you found is licensed software created by pPreetham Manjunatha, so if you need functionality that goes beyond what is currently implemented—such as computing the areas of each bounded Voronoi cell—the best approach is to contact the author directly through his GitHub page or MATLAB FileExchange profile.

Regarding your question about area computation:

Preetham’s function does not compute Voronoi regions (closed polygons). It only clips the Voronoi edge segments to a rectangular or square boundary and returns their coordinates as ‘vxClip’ and ‘vyClip’. Because the output consists of edges rather than complete cell polygons, there is no built-in method in that code to calculate the area of each region. To get areas, you would first need to reconstruct each cell’s polygon from the clipped edges, which the function itself does not provide.

As for the difference between his implementation and the other code: both perform boundary clipping of Voronoi edges, but neither produces the full polygonal regions required for area calculation. The main differences are in the internal approach to clipping and intersection handling rather than in capabilities related to region geometry.

If precise per-cell areas are essential for your application, it may help to reach out to Preetham for clarification on whether polygon reconstruction is supported or if it is something he plans to include.

Hope this helps!

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.