Floor Plan Generation with Voronoi Diagram
Introduction 馃敼
This project
reviews
and implements the paper
Free-form Floor Plan Design using Differentiable Voronoi Diagram.
In deep learning and gradient-based optimization methods, gradients are usually computed using tensor operations, but this approach is not intuitive for geometric operations.
Therefore, this project integrates tensor operations and geometric operations using
Pytorch and
Shapely.
The main difference between the paper and this project is the use of autograd.
In the paper, the
Differentiable Voronoi Diagram is used to propagate gradients,
while this project uses
Numerical Differentiation to approximate the gradients directly.
Floor plan generation with voronoi diagram
The following section describes numerical differentiation.
Numerical Differentiation
Numerical differentiation
is a method used to approximate derivatives using finite perturbation differences.
Unlike automatic differentiation used in the differentiable Voronoi Diagram in the original paper,
this approach computes derivatives by evaluating the function at nearby points.
There are three main methods for numerical differentiation, and in this project the central difference method is used to compute the gradient.
Basic methods for the numerical differentiation
Central difference method:
\[
\begin{align*}
\,\\
f'(x) &= \lim_{h \, \rightarrow \, 0} \, \frac{1}{2} \cdot \left( \frac{f(x + h) - f(h)}{h} - \frac{f(x - h) - f(h)}{h} \right)
\\\,\\
&= \lim_{h \, \rightarrow \, 0} \, \frac{1}{2} \cdot \frac{f(x + h) - f(x - h)}{h}
\\\,\\
&= \lim_{h \, \rightarrow \, 0} \, \frac{f(x + h) - f(x - h)}{2h}
\,\\
\end{align*}
\]
The value \(h \,(\text{or } dx)\) in the expression is a perturbation value that controls the accuracy of the approximation.
As \(h\) approaches zero, the numerical approximation becomes closer to the true derivative.
However, in practice, an infinitely small value cannot be used due to computational limitations and floating-point precision.
The choice of step size is important: values that are too large lead to poor approximations, while values that are too small cause numerical instability due to rounding errors.
A stable perturbation value typically ranges from \(h = 10^{-4}\) to \(h = 10^{-6}\), and in this implementation \(h = 10^{-6}\) is used.
Expression of Loss functions
In the paper, the key loss functions used to optimize floor plans consist of four main components.
The contents below are excerpted from the paper:
-
Wall loss:
As the unconstrained Voronoi diagram typically produce
undesirable fluctuations in the wall orientations,
we design a tailored loss function to regularize the wall complexity.
Inspired by the Cubic Stylization,
we regularize the \(\mathcal{L}_1\) norm of the wall length.
\(L_1\) norm is defined as \(v_x + v_y\) (norm of \(x\) + norm of \(y\)),
therefore the \(\mathcal{L}_{\text{wall}}\) has the minimal when vector \(\mathbb{v}_j - \mathbb{v}_i\) is vertical or horizontal.
\[
\,\\
\mathcal{L}_{\text{wall}} = w_{\text{wall}} \sum_{(v_i, v_j) \, \in \, \mathcal{E}} ||\, \mathbb{v}_i - \mathbb{v}_j \,||_{L1}
\,\\
\]
where \(\mathcal{E}\) denotes the set of edges of the Voronoi cells between two
adjacent rooms and the \(\mathbb{v}_i\) and \(\mathbb{v}_j\) denote the Voronoi vertices belonging to the edge.
-
Area loss: The area of each room is specified by the user.
We minimize the quadratic difference between the current room areas and
the user-specified targets. Here, \(\bar{A}_r\) denotes the target area for the room \(r\).
\[
\,\\
\mathcal{L}_{\text{area}} = w_{\text{area}} \sum_{r=1}^{\#Room} ||\, A_r(\mathcal{V}) - \bar{A}_r \,||^2
\,\\
\]
-
Lloyd loss: To regulate the site density, we design a loss function inspired by the Lloyd's algorithm.
Here, \(\mathbb{c}_i \) denotes the centroid of the \(i\)-th Voronoi cell.
This is useful for attracting these exterior sites inside \(\Omega\).
\[
\,\\
\mathcal{L}_{\text{Lloyd}} = w_{\text{Lloyd}} \sum_{i=1}^N ||\, \mathbb{s}_i - \mathbb{c}_i \,||^2
\,\\
\]
-
Topology loss: We design the topology loss such that each room is
a single connected region, and the specified connections between
rooms are achieved. We move the site to satisfy the desired topology by setting the goal position
\(\mathbb{t}_i\) for each site
\(\mathbb{s}_i\) as
\[
\,\\
\mathcal{L}_{\text{topo}} = w_{\text{topo}} \sum_{i=1}^N ||\, \mathbb{s}_i - \mathbb{t}_i \,||^2
\,\\
\]
The goal position
\(\mathbb{t}_i\)
can be automatically computed as the nearest
site to the site from the same group.
For each room, we
first group the sites belonging to that room into groups of adjacent
sites. If multiple groups are present, that is, a room is split into
separated regions,
we set the target position of the site \(\mathbb{t}_i\)
as the
nearest site to that group.
Implementation of loss functions
As mentioned in the Introduction, the loss functions above are implemented for forward propagation using Shapely and Pytorch, as shown below.
The total
loss is defined as a weighted sum of the individual losses, and this loss is used to update the Voronoi sites and generate the floor plan.
\[
\,\\
\begin{align*}
\mathcal{S}^{*} &= \arg \min_{\mathcal{S}} \mathcal{L}(\mathcal{S}, \mathcal{V}(\mathcal{S})) \\
\mathcal{L} &= \mathcal{L}_{\text{wall}} + \mathcal{L}_{\text{area}} + \mathcal{L}_{\text{fix}} + \mathcal{L}_{\text{topo}} + \mathcal{L}_{\text{Lloyd}}
\end{align*}
\,\\
\]
class FloorPlanLoss(torch.autograd.Function):
@staticmethod
def compute_wall_loss(rooms_group: List[List[geometry.Polygon]], w_wall: float = 1.0):
loss_wall = 0.0
for room_group in rooms_group:
room_union = ops.unary_union(room_group)
if isinstance(room_union, geometry.MultiPolygon):
room_union = list(room_union.geoms)
else:
room_union = [room_union]
for room in room_union:
t1 = torch.tensor(room.exterior.coords[:-1])
t2 = torch.roll(t1, shifts=-1, dims=0)
loss_wall += torch.abs(t1 - t2).sum().item()
for interior in room.interiors:
t1 = torch.tensor(interior.coords[:-1])
t2 = torch.roll(t1, shifts=-1, dims=0)
loss_wall += torch.abs(t1 - t2).sum().item()
loss_wall = torch.tensor(loss_wall)
loss_wall *= w_wall
return loss_wall
@staticmethod
def compute_area_loss(
cells: List[geometry.Polygon],
target_areas: List[float],
room_indices: List[int],
w_area: float = 1.0,
):
current_areas = [0.0] * len(target_areas)
for cell, room_index in zip(cells, room_indices):
current_areas[room_index] += cell.area
current_areas = torch.tensor(current_areas)
target_areas = torch.tensor(target_areas)
area_difference = torch.abs(current_areas - target_areas)
loss_area = torch.sum(area_difference)
loss_area **= 2
loss_area *= w_area
return loss_area
@staticmethod
def compute_lloyd_loss(cells: List[geometry.Polygon], sites: torch.Tensor, w_lloyd: float = 1.0):
valids = [(site.tolist(), cell) for site, cell in zip(sites, cells) if not cell.is_empty]
valid_centroids = torch.tensor([cell.centroid.coords[0] for _, cell in valids])
valid_sites = torch.tensor([site for site, _ in valids])
loss_lloyd = torch.norm(valid_centroids - valid_sites, dim=1).sum()
loss_lloyd **= 2
loss_lloyd *= w_lloyd
return loss_lloyd
@staticmethod
def compute_topology_loss(rooms_group: List[List[geometry.Polygon]], w_topo: float = 1.0):
loss_topo = 0.0
for room_group in rooms_group:
room_union = ops.unary_union(room_group)
if isinstance(room_union, geometry.MultiPolygon):
largest_room, *_ = sorted(room_union.geoms, key=lambda r: r.area, reverse=True)
loss_topo += len(room_union.geoms)
for room in room_group:
if not room.intersects(largest_room) and not room.is_empty:
loss_topo += largest_room.centroid.distance(room)
loss_topo = torch.tensor(loss_topo)
loss_topo **= 2
loss_topo *= w_topo
return loss_topo
( ... )
@staticmethod
def forward(
ctx: FunctionCtx,
sites: torch.Tensor,
boundary_polygon: geometry.Polygon,
target_areas: List[float],
room_indices: List[int],
w_wall: float,
w_area: float,
w_lloyd: float,
w_topo: float,
w_bb: float,
w_cell: float,
save: bool = True,
) -> torch.Tensor:
cells = []
walls = []
sites_multipoint = geometry.MultiPoint([tuple(point) for point in sites.detach().numpy()])
raw_cells = list(shapely.voronoi_polygons(sites_multipoint, extend_to=boundary_polygon).geoms)
for cell in raw_cells:
intersected_cell = cell.intersection(boundary_polygon)
intersected_cell_iter = [intersected_cell]
if isinstance(intersected_cell, geometry.MultiPolygon):
intersected_cell_iter = list(intersected_cell.geoms)
for intersected_cell in intersected_cell_iter:
exterior_coords = torch.tensor(intersected_cell.exterior.coords[:-1])
exterior_coords_shifted = torch.roll(exterior_coords, shifts=-1, dims=0)
walls.extend((exterior_coords - exterior_coords_shifted).tolist())
cells.append(intersected_cell)
cells_sorted = []
raw_cells_sorted = []
for site_point in sites_multipoint.geoms:
for ci, (cell, raw_cell) in enumerate(zip(cells, raw_cells)):
if raw_cell.contains(site_point):
cells_sorted.append(cell)
cells.pop(ci)
raw_cells_sorted.append(raw_cell)
raw_cells.pop(ci)
break
rooms_group = [[] for _ in torch.tensor(room_indices).unique()]
for cell, room_index in zip(cells_sorted, room_indices):
rooms_group[room_index].append(cell)
loss_wall = torch.tensor(0.0)
if w_wall > 0:
loss_wall = FloorPlanLoss.compute_wall_loss(rooms_group, w_wall=w_wall)
loss_area = torch.tensor(0.0)
if w_area > 0:
loss_area = FloorPlanLoss.compute_area_loss(cells_sorted, target_areas, room_indices, w_area=w_area)
loss_lloyd = torch.tensor(0.0)
if w_lloyd > 0:
loss_lloyd = FloorPlanLoss.compute_lloyd_loss(cells_sorted, sites, w_lloyd=w_lloyd)
loss_topo = torch.tensor(0.0)
if w_topo > 0:
loss_topo = FloorPlanLoss.compute_topology_loss(rooms_group, w_topo=w_topo)
loss_bb = torch.tensor(0.0)
if w_bb > 0:
loss_bb = FloorPlanLoss.compute_bb_loss(rooms_group, w_bb=w_bb)
loss_cell_area = torch.tensor(0.0)
if w_cell > 0:
loss_cell_area = FloorPlanLoss.compute_cell_area_loss(cells_sorted, w_cell=w_cell)
if save:
ctx.save_for_backward(sites)
ctx.room_indices = room_indices
ctx.target_areas = target_areas
ctx.boundary_polygon = boundary_polygon
ctx.w_wall = w_wall
ctx.w_area = w_area
ctx.w_lloyd = w_lloyd
ctx.w_topo = w_topo
ctx.w_bb = w_bb
ctx.w_cell = w_cell
loss = loss_wall + loss_area + loss_lloyd + loss_topo + loss_bb + loss_cell_area
return loss, [loss_wall, loss_area, loss_lloyd, loss_topo, loss_bb, loss_cell_area]
Since the loss functions were implemented using Shapely in Python, there are some differences compared to the original implementation.
Backward with numerical differentiation
Numerical differentiation is not efficient in terms of computational performance, because it requires multiple function evaluations at nearby points to approximate derivatives.
As shown in the
backward method below, computational performance depends on the number of sites.
Therefore, Python's built-in multiprocessing module was used to improve the performance of
backward propagation.
@staticmethod
def _backward_one(args):
(
sites,
i,
j,
epsilon,
boundary_polygon,
target_areas,
room_indices,
w_wall,
w_area,
w_lloyd,
w_topo,
w_bb,
w_cell,
) = args
perturbed_sites_pos = sites.clone()
perturbed_sites_neg = sites.clone()
perturbed_sites_pos[i, j] += epsilon
perturbed_sites_neg[i, j] -= epsilon
loss_pos, _ = FloorPlanLoss.forward(
None,
perturbed_sites_pos,
boundary_polygon,
target_areas,
room_indices,
w_wall,
w_area,
w_lloyd,
w_topo,
w_bb,
w_cell,
save=False,
)
loss_neg, _ = FloorPlanLoss.forward(
None,
perturbed_sites_neg,
boundary_polygon,
target_areas,
room_indices,
w_wall,
w_area,
w_lloyd,
w_topo,
w_bb,
w_cell,
save=False,
)
return i, j, (loss_pos - loss_neg) / (2 * epsilon)
@runtime_calculator
@staticmethod
def backward(ctx: FunctionCtx, _: torch.Tensor, __):
sites = ctx.saved_tensors[0]
room_indices = ctx.room_indices
target_areas = ctx.target_areas
boundary_polygon = ctx.boundary_polygon
w_wall = ctx.w_wall
w_area = ctx.w_area
w_lloyd = ctx.w_lloyd
w_topo = ctx.w_topo
w_bb = ctx.w_bb
w_cell = ctx.w_cell
epsilon = 1e-6
grads = torch.zeros_like(sites)
multiprocessing_args = [
(
sites,
i,
j,
epsilon,
boundary_polygon,
target_areas,
room_indices,
w_wall,
w_area,
w_lloyd,
w_topo,
w_bb,
w_cell,
)
for i in range(sites.size(0))
for j in range(sites.size(1))
]
with multiprocessing.Pool(processes=multiprocessing.cpu_count()) as pool:
results = pool.map(FloorPlanLoss._backward_one, multiprocessing_args)
for i, j, grad in results:
grads[i, j] = grad
return grads, None, None, None, None, None, None, None, None, None, None
Initializing parameters
In optimization problems, the initial parameters significantly affect the final results.
Firstly, I initialized the Voronoi diagram's
sites such that the sites were generated at the center of a given floor plan:
-
Random Sites Generation: Generate initial random sites using uniform distribution.
-
Moving to Center of Boundary: Shift all sites so they are centered within the floor plan boundary.
-
Outside Sites Adjustment: Adjust any sites that fall outside the boundary by moving them inward.
-
Voronoi Diagram: Generate Voronoi diagram using sites.
Process of parameters initialization
Secondly, I used the
KMeans clustering algorithm to assign cell indices to each site.
The distance-based KMeans algorithm groups sites by spatial proximity, which helps ensure that rooms are formed from adjacent cells.
By clustering the sites in advance, initial room assignments become spatially coherent,
which reduces the possibility of disconnected room regions during optimization.
Using this approach, the optimizer converges more stably. An example is shown below.
Floor plan generation on 300 iterations
From the left, optimization without KMeans 路 optimization with KMeans
As shown in the figure above, KMeans makes the loss change more smoothly and converge faster.
Without KMeans, the optimization process shows unstable behavior and often produces disconnected rooms.
In contrast, when KMeans is used for initial room assignments, the optimization maintains spatial coherence throughout the process, leading to:
-
faster convergence to the target room areas,
-
more stable wall alignments, and
-
a reduced possibility of rooms splitting into disconnected regions.
This improvement in optimization stability is particularly important for complex floor plans with multiple rooms and specific area requirements.
Experiments
Future works
-
Set entrance: In the paper, to set entrance of the plan it uses \(\mathcal{L}_{\text{fix}}\) loss function.
-
Graph-based contraint: In the paper, to set and ensure the rooms' adjacencies it uses graph-based constraint.
-
Improve computational performance: Optimize the code to run faster, for example by porting it to another language or by implementing a differentiable Voronoi diagram.
-
Handle deadspaces: Set the loss function for deadspace \(\mathcal{L}_{\text{deadspace}}\) to exclude infeasible plans.
-
Following boundary axis: Align walls to follow the axis of a given boundary instead of global X, Y (Replacing \(\mathcal{L}_{\text{wall}}\)).