The greedy approach wouldn't be that hard to implement, and doesn't require polygonal approximations.

Just think of $P$ as the space between zero-curvature disks. Then at each step of the algorithm, you have a region bounded by disks, and you want to find the largest disk that is tangent to three boundary disks, and not intersecting any others. After this, you are left with three (generically, possibly more if things are degenerate) regions to add to your stack.

You want to this depth-first, to avoid using too much memory, so you should stop going down a branch if the remaining uncovered space is $<3^{-d} \epsilon$, where $d$ is the depth.

The greedy approach wouldn't be that hard to implement, and doesn't require polygonal approximations.

Just think of $P$ as the space between zero-curvature disks. Then at each step of the algorithm, you have a region bounded by disks, and you want to find the largest disk that is tangent to three boundary disks, and not intersecting any others. After this, you are left with three (generically, possibly more if things are degenerate) regions to add to your stack.

You want to this depth-first, to avoid using too much memory, so you should stop going down a branch if the region has area $<3^{-d} \epsilon$, where $d$ is the depth. Alternatively, stop if the region can be circumscribed by a disk such that the overlap created is $<3^{-d}\delta$.