diff --git a/cpp/locality/Voronoi.cc b/cpp/locality/Voronoi.cc
index 35feb907f..35e58c733 100644
--- a/cpp/locality/Voronoi.cc
+++ b/cpp/locality/Voronoi.cc
@@ -16,7 +16,7 @@
 namespace freud { namespace locality {
 
 // Voronoi calculations should be kept in double precision.
-void Voronoi::compute(const freud::locality::NeighborQuery* nq)
+void Voronoi::compute(const freud::locality::NeighborQuery* nq, const double* radii)
 {
     const auto box = nq->getBox();
     const auto n_points = nq->getNPoints();
@@ -39,13 +39,18 @@ void Voronoi::compute(const freud::locality::NeighborQuery* nq)
     const int voro_blocks_y = int(box.getLy() * block_scale + 1);
     const int voro_blocks_z = int(box.getLz() * block_scale + 1);
 
-    voro::container_periodic container(v1.x, v2.x, v2.y, v3.x, v3.y, v3.z, voro_blocks_x, voro_blocks_y,
-                                       voro_blocks_z, 3);
+    voro::container_periodic_poly container(v1.x, v2.x, v2.y, v3.x, v3.y, v3.z, voro_blocks_x, voro_blocks_y,
+                                            voro_blocks_z, 3);
 
     for (size_t query_point_id = 0; query_point_id < n_points; query_point_id++)
     {
+        // If an array of radii is provided, those radii are used to compute
+        // the power diagram (also called the radical Voronoi tessellation). If
+        // a nullptr is provided, then the points are assumed to have zero
+        // radius.
         vec3<double> query_point((*nq)[query_point_id]);
-        container.put(query_point_id, query_point.x, query_point.y, query_point.z);
+        double radius = (radii != nullptr) ? radii[query_point_id] : 0.0;
+        container.put(query_point_id, query_point.x, query_point.y, query_point.z, radius);
     }
 
     voro::voronoicell_neighbor cell;
diff --git a/cpp/locality/Voronoi.h b/cpp/locality/Voronoi.h
index 01ecad957..6bf69307b 100644
--- a/cpp/locality/Voronoi.h
+++ b/cpp/locality/Voronoi.h
@@ -19,7 +19,7 @@ class Voronoi
     // default constructor
     Voronoi() : m_neighbor_list(std::make_shared<NeighborList>()) {}
 
-    void compute(const freud::locality::NeighborQuery* nq);
+    void compute(const freud::locality::NeighborQuery* nq, const double* radii = nullptr);
 
     std::shared_ptr<NeighborList> getNeighborList() const
     {
diff --git a/freud/_locality.pxd b/freud/_locality.pxd
index 8b5518b90..128397d55 100644
--- a/freud/_locality.pxd
+++ b/freud/_locality.pxd
@@ -137,7 +137,8 @@ cdef extern from "PeriodicBuffer.h" namespace "freud::locality":
 cdef extern from "Voronoi.h" namespace "freud::locality":
     cdef cppclass Voronoi:
         Voronoi()
-        void compute(const NeighborQuery*) nogil except +
+        void compute(const NeighborQuery*,
+                     const double*) nogil except +
         vector[vector[vec3[double]]] getPolytopes() const
         const freud.util.ManagedArray[double] &getVolumes() const
         shared_ptr[NeighborList] getNeighborList() const
diff --git a/freud/locality.pyx b/freud/locality.pyx
index 29d0b622c..b5ba76d3e 100644
--- a/freud/locality.pyx
+++ b/freud/locality.pyx
@@ -1162,17 +1162,28 @@ cdef class Voronoi(_Compute):
     def __dealloc__(self):
         del self.thisptr
 
-    def compute(self, system):
+    def compute(self, system, radii=None):
         r"""Compute Voronoi diagram.
 
         Args:
             system:
                 Any object that is a valid argument to
                 :class:`freud.locality.NeighborQuery.from_system`.
+            radii ((:math:`N_{points}`) :class:`numpy.ndarray`):
+                An array of radii for each point in the system. If provided,
+                the power diagram (also called the radical Voronoi
+                tessellation) will be computed (Default value = :code:`None`,
+                which gives the Voronoi diagram).
         """
         cdef NeighborQuery nq = NeighborQuery.from_system(system)
-        self.thisptr.compute(nq.get_ptr())
         self._box = nq.box
+        cdef double* l_radii_ptr = NULL
+        cdef double[::1] l_radii
+        if radii is not None:
+            l_radii = freud.util._convert_array(
+                radii, shape=(len(nq.points),), dtype=np.float64)
+            l_radii_ptr = &l_radii[0]
+        self.thisptr.compute(nq.get_ptr(), l_radii_ptr)
         return self
 
     @_Compute._computed_property