reading boundaries from *_physical_region.xml file

When dealing with custom meshes with subdomains, my typical workflow is to generate a mesh in Gmsh, convert it to XML format using dolfin-convert, and then read the mesh from the two files mymesh.xml and mymesh_physical_region.xml. The second file gives me the subdomains, e.g.,

mesh = Mesh('mymesh.xml')
subdomains = MeshFunction('size_t', mesh, 'mymesh_physical_region.xml')

From there, I typically go ahead and define the volumes,

dx = Measure('dx')[subdomains]

What I'd like to do now I also define subdomain boundary "volumes", e.g.,

ds = Measure('ds')[boundaries]

How can I generate the corresponding `boundaries`?