SHAPE AWARE AUTOMATIC REGION-OF-INTEREST SUBDIVISIONS

The fast marching algorithm [9] was then utilized to locate one of the object’s ends (Fig. 2 upper right), and then a second time (Fig. 2 lower left) to aid in computing a discrete, shortest path curve that defines the region’s centerline (Fig. 2 lower right). Briefly, the fast marching algorithm begins at a source (in this case a single voxel) and simulates a wave emanating in all directions. The propagating wave is controlled by a weighting map (i.e., the distance map) which tells the algorithm how fast the wave travels through each voxel. As the wave passes over each voxel, “how long” it took the wave to get to that voxel is recorded, termed the “arrival time”. Therefore, the arrival time is a function of both the distance a particular voxel is from the source of the wave front, as well as the weighting function that the wave front experiences along the way.

In $\mathbb{R}^{2}$ we begin with a potential function $V(x)$ that weights a region’s voxels according to their distance from the region’s border. The weighted geodesic distance between two points $x_{0},x_{1}\in\mathbb{R}^{2}$ is given by[10]

Shortest paths can then be extracted from any voxel to the wave front’s source voxel (thus finding the geodesic curve $\gamma$ from $x_{0}$ to $x_{1}$)[12]. A discrete gradient descent method was utilized to compute this path[13].

Fig. 2 depicts snapshots along the centerline defining algorithm. The distance map of the region is first computed (Fig. 2 upper left) and the maximum value found is used as the starting point for the first fast marching wave (Fig. 2 upper right). The maximum value found from this fast marching wave is then used to propagate the second and final fast marching wave (which is weighted by the distance map to the sixth power so that the wave travels fastest in the region’s center), as shown in the lower left panel of Fig. 2. Finally a discrete gradient descent is performed starting from the maximum voxel in the image created by the second fast marching wave. This shortest path defines the object’s centerline as shown in the lower right panel of Fig. 2 overlaid on the original binary mask.

To divide the ROI into subdivisions, cross-sectional areas, or line-profiles were calculated at each centerline voxel. In our current example we illustrate the division of a single ROI into 16 subregions. First, 15 evenly spaced voxels are labeled in the ordered list of centerline voxels. For each of the 15 voxels, the normal vector, or direction of the centerline, is computed from the two adjacent voxels. The adjacent voxels of ($x_{i}$,$y_{i}$) are labelled ($x_{i-1}$,$y_{i-1}$) and ($x_{i+1}$,$y_{i+1}$). The normal vector is then calculated as

Three example cross-sections, or line profiles cutting through the region, are shown in the first three panels of Fig. 3. Separating the object by these cuts, the individual regions can be delineated by a connected component analysis. To label the individual regions, the algorithm iterated through each cut, and added the removed region to the existing image. The result was a mask with regions labeled from 1 to 16, thereby defining the SROIs (right panel of Fig. 3).

In the case of certain types of analysis, it may also be desirable to keep the area (or volume) of each subregion the same. To create equal area regions an additional step was implemented. First, the largest region was considered (call $R_{l}$). Next, the two neighboring regions’ areas were computed and the smaller of the two was chosen (call $R_{s}$). All voxels of $R_{l}$ bordering with $R_{s}$ were changed to be labeled as part of $R_{s}$. This process was repeated (until no new $R_{l}$’s were chosen) in order to perform a rough area correction that retains the shape of the individual regions (only changing the width). Finally the small differences between areas were corrected by exchanging bordering voxels from larger regions with those of neighboring smaller regions. Since the area of the ROI will most likely not be a multiple of the number of SROIs to be generated, a final automation step removed a few outer voxels of the last region (starting from those with a maximum value in the second fast marching wave image) so that each of the regions have equal areas.