diff --git a/pycoloquinte/module.cpp b/pycoloquinte/module.cpp
index 4a7c8a6..e4529e6 100644
--- a/pycoloquinte/module.cpp
+++ b/pycoloquinte/module.cpp
@@ -318,6 +318,10 @@ Construct a circuit.
            py::arg("expansion_factor"), py::arg("max_density") = 1.0,
            py::arg("row_side_margin") = 0.0,
            "Expand the standard cells by an individual factor")
+      .def("compute_cell_expansion", &Circuit::computeCellExpansion,
+           py::arg("congestion_map"), py::arg("fixed_penalty") = 0.0,
+           py::arg("penalty_factor") = 1.0,
+           "Compute an expansion factor given a congestion map")
       .def("check", &Circuit::check, "Check the datastructure")
       .def("report", &Circuit::report)
       .def("export_ispd", &Circuit::exportIspd,
diff --git a/src/coloquinte.cpp b/src/coloquinte.cpp
index 3e11d8c..8e5302c 100644
--- a/src/coloquinte.cpp
+++ b/src/coloquinte.cpp
@@ -727,14 +727,28 @@ void Circuit::expandCellsByFactor(const std::vector<float> &expansionFactor,
 }
 
 std::vector<float> Circuit::computeCellExpansion(
-    const std::vector<std::pair<Rectangle, float> > &congestionMap) const {
+    const std::vector<std::pair<Rectangle, float> > &congestionMap,
+    float fixedPenalty, float penaltyFactor) const {
+  if (fixedPenalty < 0.0f || penaltyFactor < 1.0f) {
+    throw std::runtime_error(
+        "For cell expansion, fixed penalty should be non-negative, and penalty "
+        "factor should be at least 1");
+  }
   std::vector<std::pair<Rectangle, float> > expansionMap;
   for (auto [r, c] : congestionMap) {
-    if (c > 0.0f) {
-      expansionMap.emplace_back(r, c + 1.0f);
+    if (c > 1.0f) {
+      expansionMap.emplace_back(
+          r, (c - 1.0f) * penaltyFactor + fixedPenalty + 1.0);
     }
   }
-  //std::sort(expansionMap.begin(), expansionMap.end());
+
+  // Sort to allow binary search later
+  std::sort(
+      expansionMap.begin(), expansionMap.end(),
+      [](const CongestionRegion &a, const CongestionRegion &b) -> bool {
+        return a.first.minX < b.first.minX ||
+               (a.first.minX == b.first.minX && a.first.minY < b.first.minY);
+      });
 
   // Now analyze the expansion for each cell; use the maximum of the expansion
   // maps it intersects
@@ -745,6 +759,7 @@ std::vector<float> Circuit::computeCellExpansion(
     } else {
       Rectangle place = placement(i);
       float expansion = 1.0;
+      // TODO: use binary search here instead of dumb iteration
       for (auto [r, e] : expansionMap) {
         if (r.intersects(place)) {
           expansion = std::max(expansion, e);
diff --git a/src/coloquinte.hpp b/src/coloquinte.hpp
index 12f2da8..1eb9305 100644
--- a/src/coloquinte.hpp
+++ b/src/coloquinte.hpp
@@ -1024,15 +1024,24 @@ class Circuit {
    * coordinates from the original circuit.
    */
 
+  /**
+   * @brief Representation of a congested region, as a congestion factor
+   */
+  using CongestionRegion = std::pair<Rectangle, float>;
+
   /**
    * @brief Compute the cell expansion required for a given congestion map
    *
    * @param congestionMap Region and associated congestion, as a proportion of
-   * the available routing resources (0.0 -> not congested; 0.5 -> 50% over
-   * routing capacity).
+   * the available routing capacity:
+   *   * 0.0 -> empty
+   *   * 0.5 -> 50% occupancy
+   *   * 1.0 -> 100% occupancy
+   *   * 1.1 -> 110% occupancy (CONGESTED)
    */
   std::vector<float> computeCellExpansion(
-      const std::vector<std::pair<Rectangle, float> > &congestionMap) const;
+      const std::vector<CongestionRegion> &congestionMap,
+      float fixedPenalty = 0.0, float penaltyFactor = 1.0) const;
 
   /**
    * @brief Return a brief description of the circuit
diff --git a/test/test_expansion.cpp b/test/test_expansion.cpp
index 45a1f11..9af0975 100644
--- a/test/test_expansion.cpp
+++ b/test/test_expansion.cpp
@@ -62,13 +62,18 @@ BOOST_AUTO_TEST_CASE(TestFactorExpansionMultiple) {
 }
 
 BOOST_AUTO_TEST_CASE(TestExpansionComputation) {
-  Circuit circuit(1);
-  circuit.setCellWidth({20});
-  circuit.setCellHeight({10});
-  circuit.setCellX({40});
-  circuit.setCellY({-100});
+  Circuit circuit(4);
+  circuit.setCellWidth({20, 20, 20, 20});
+  circuit.setCellHeight({10, 10, 10, 10});
+  circuit.setCellX({40, 9, 10, -5});
+  circuit.setCellY({-100, 19, 20, 19});
   std::vector<std::pair<Rectangle, float> > congestionMap;
+  congestionMap.emplace_back(Rectangle({0, 10, 10, 20}), 1.2);
+  congestionMap.emplace_back(Rectangle({-5, 0, 10, 20}), 1.5);
   std::vector<float> expansion = circuit.computeCellExpansion(congestionMap);
   BOOST_CHECK_EQUAL(expansion.size(), circuit.nbCells());
   BOOST_CHECK_EQUAL(expansion[0], 1.0f);
+  BOOST_CHECK_EQUAL(expansion[1], 1.2f);
+  BOOST_CHECK_EQUAL(expansion[2], 1.0f);
+  BOOST_CHECK_EQUAL(expansion[3], 1.5f);
 }