1. Algorithm description
The developed algorithm herein is scripted in Python and executed using Python version 3.8.10. The primary libraries necessary to construct the algorithm are Pydicom, NumPy, OpenCV and SciPy. Pydicom library enables the reading of DICOM (Digital Imaging and Communications in Medicine) files that hold dose distributions [8], which will be manipulated using NumPy library. The OpenCV library is employed for various stages of image registration. In order to compare the resulting distributions, the linear least-square regression method was employed through the SciPy library.
The algorithm applies the image registration approach for dose distributions comparison using ten functions alongside the main program. These functions detect keypoints from the isodose contours of dose distributions being compared and compute their corresponding descriptors [6–8]. Then, they match keypoints of both distributions by computing the distance between their descriptors and aligning keypoints with the minimum distance [6][9]. Finally, they estimate an affine transformation matrix that transforms one dose distribution to the coordinates of the other, as depicted in Fig. 1 [9–11].
Within the main program, the algorithm reads two DICOM files containing the dose distributions. It labels the smaller distribution as the evaluated distribution and the larger one as the reference distribution. The transformation matrix is applied to the evaluated dose distribution, resulting in a new distribution with the same dimensions as the reference distribution (Fig. 2). Consequently, dose comparison can be accomplished by examining the slope between the transformed distribution and the reference distribution.
1.1. Keypoints detection
Keypoints are generated based on convexity defects in isodose contours of each dose distribution. Convexity defects of a contour are deviations of the contour from its convex hull, which is the minimum convex boundary that may envelope the contour. These defects are produced by detecting isodose contours for each dose distribution, creating convex hull for each contour and extracting points that define each convexity defect (the starting point, the ending point and the farthest point from the hull) [12–13].
Contour detection requires pre-processing of dose distributions that involves thresholding and applying the Laplacian filter [7–8]. The resulting binary matrix (Fig. 3) is used as an argument in the function ‘findContours’ of OpenCV library to find each isodose contour points [14].
Before creating the contours’ convex hulls, each isodose contour of the evaluated dose distribution is compared with its corresponding contour in the reference dose distribution, and dissimilar contours are removed. The 'convexHull' and 'convexityDefects' functions within the OpenCV library are utilized to identify convex hulls and subsequently extract the points corresponding to their convexity defects. The three points defining each convexity defect are considered as a single keypoint (Fig. 4)[14].
1.2. Keypoints description
Each keypoint's descriptor is a vector with three components. The first component (D) is the Euclidean distances between the contour's centroid and the three points of the keypoint, divided by the minimum enclosing circle's diameter of the contour. This component includes also the Euclidean distances between the three points, each divided by the minimum enclosing circle's radius of the contour. The second component is the histogram of oriented gradients (HOG) of the convexity defect, and the third is the curvature histogram (K) of the convexity defect curve.
The histogram of oriented gradients component is computed using gradient orientations (Eq. 6) and gradient magnitudes (Eq. 5) of convexity defect points. Histogram bins represent gradient orientations, and bin values represent the corresponding weighted gradient magnitudes [15]. These gradients are estimated using Scharr operator as shown in the equations 1, 2, 3 and 4.
\({S}_{x}=\left[\begin{array}{ccc}3& 0& -3\\ 10& 0& -10\\ 3& 0& -3\end{array}\right]\)
|
(1)
|
\({S}_{y}=\left[\begin{array}{ccc}3& 10& 3\\ 0& 0& 0\\ -3& -10& -3\end{array}\right]\)
|
(2)
|
\({G}_{x}=D*{S}_{x}\)
|
(3)
|
\({ G}_{y}=D*{S}_{y}\)
|
(4)
|
\(G=\sqrt{{{G}_{x}}^{2}+{{G}_{y}}^{2}}\)
|
(5)
|
\(Ɵ=atan\left(\frac{{G}_{y}}{{G}_{x}}\right)\)
|
(6)
|
Where \({S}_{x}\) and \({S}_{y}\) are the horizontal and vertical operators of Scharr operator, respectively; while \(G\) is the gradient magnitude and \(Ɵ\) is its orientation [16].
After dividing the range of gradient orientations into a fixed number of bins (40 by default), the contribution weight of each gradient magnitude to the histogram is calculated based on the proximity of the gradient's orientation to the corresponding bin center. For example, if the orientation of a point was 26 for a bin size of 10, it contributes to the bin [20–30] with 0.4 and to the bin [30–40] with 0.6. Weighted gradient magnitudes of each bin are, then, accumulated to build the histogram of oriented gradients.
The curvature histogram is generated by accumulating the weighted frequencies of curvature values into their corresponding bins in the histogram. The curvature \(K\) of each point in the convexity defect is computed using the Eq. 7.
$$\varvec{K}=\frac{\left|{\varvec{y}}^{\varvec{{\prime }}}{\varvec{x}}^{\varvec{{\prime }}\varvec{{\prime }}}-\varvec{x}\varvec{{\prime }}\varvec{y}\varvec{{\prime }}\varvec{{\prime }}\right|}{{\left({\varvec{x}\varvec{{\prime }}}^{2}+{\varvec{y}\varvec{{\prime }}}^{2}\right)}^{3/2}}$$
7
Where \(x\) and \(y\) represent the coordinates of the defect point along the x-axis and the y-axis, respectively; \({x}^{{\prime }}\) and \({y}^{{\prime }}\) denote their first derivatives, while \({x}^{{\prime }{\prime }}\) and \({y}^{{\prime }{\prime }}\)represent their second derivatives [7].
Each descriptor's component is normalized by its elements' sum, and the total descriptor is normalized by the normalized components' sum, as described in the equations 8, 9, 10 and 11.
\(N\_{D}_{i}= \frac{{D}_{i}}{\sum _{i}{D}_{i}}\)
|
(8)
|
\({N\_HOG}_{i}= \frac{{HOG}_{i}}{\sum _{i}{HOG}_{i}}\)
|
(9)
|
\({K}_{i}= \frac{{K}_{i}}{\sum _{i}{K}_{i}}\)
|
(10)
|
\({N\_Descp}_{i}= \frac{{Descp}_{i}}{\sum _{i}{D}_{i}+\sum _{i}{HOG}_{i}+\sum _{i}{K}_{i}}\)
|
(11)
|
Where \({D}_{i}\), \({HOG}_{i}\) and \({K}_{i}\) are the distance component, the oriented gradients component and the curvature component for each keypoint, respectively; and \({Descp}_{i}\) is the vector containing these components. \(N\_{D}_{i}\), \({N\_HOG}_{i}\), \({K}_{i}\) and \(N\_Descp\) are their normalized counterparts.
1.3. Keypoints matching
During the matching step, the keypoints of each isodose contour in the evaluated dose distribution are aligned with their counterparts in the reference dose distribution (Fig. 5). This alignment is achieved by comparing each keypoint in the contour with the smallest number of keypoints to a set of keypoints in the counterpart contour of the other distribution. The number of keypoints in the second contour that will be compared with the keypoints of the first contour is determined by the ratio k. This ratio is obtained by dividing the number of keypoints in the second contour by the number of keypoints in the first contour.
The algorithm then calculates the Manhattan distance between the descriptor of the ith keypoint in the first contour and a range from \(⌊ki⌋\) to \(⌈k(i+1)⌉\) of descriptors in the second contour. The keypoint pair that corresponds to the descriptors with the minimum distance is matched, and the matching information is stored in a \(m\times n\) matrix, where \(m\) is the number of keypoints in the first contour and \(n\) is the number of keypoints in the second contour (\(m\le n\)).
If the ith keypoint in the first contour is matched with the jth keypoint in the second contour, then the (i,j) element of the matrix is set to 1; otherwise, it is set to 0. No keypoint in the first contour can have more than one match in the second. On the other hand, keypoints in the second contour can match more than one keypoint in the first contour because they can be part of two different sets: one for comparing with a keypoint in the first contour and another for comparing with the next keypoint.
$$\left[\begin{array}{ccc}1& \cdots & 0\\ ⋮& \begin{array}{ccc}\ddots & 0& 0\\ 0& 0& 0\\ 0& 0& \ddots \end{array}& ⋮\\ 0& \cdots & 0\end{array}\right]$$a. Matching matrix where the first keypoints in both dose distributions are matched.
$$\left[\begin{array}{ccc}0& \cdots & 0\\ ⋮& \begin{array}{ccc}\ddots & 1& 0\\ 0& 1& 0\\ 0& 0& \ddots \end{array}& ⋮\\ 0& \cdots & 0\end{array}\right]$$
b. Matching matrix that violates one-to-one keypoint matching constraint.
$$\left[\begin{array}{ccc}0& \cdots & 0\\ ⋮& \begin{array}{ccc}\ddots & 1& 0\\ 1& 0& 0\\ 0& 0& \ddots \end{array}& ⋮\\ 0& \cdots & 0\end{array}\right]$$c. Matching matrix that violates order-preserving matching constraint.
Our matching approach adheres to two constraints: one-to-one keypoint matching and order-preserving matching. The first constraint guarantees that each keypoint in the first contour is matched with a single keypoint in the second contour, and vice versa. The second constraint ensures that the order of matched keypoints along the contours is maintained, so that no two matches intersect each other.
If any matches do not meet these constraints, only the match with the minimum distance is kept. The rest are set to 0 until a new match is found by adjusting the range of keypoints in the second contour to be compared with.
1.4. Applying the transformation
After matching all the isodose contours in both distributions, two lists of keypoints are obtained. Each keypoint in the first list corresponds to the keypoint in the second list with the same order. These two lists are used as arguments to the ‘estimateAffine2D’ function from OpenCV library, which estimates the affine transformation matrix between the reference distribution and the evaluated distribution. By applying this transformation to the evaluated distribution using ‘warpAffine’ function from OpenCV, a matrix that has the same coordinates as the reference distribution is acquired. This matrix can be used in the comparison step, in which the algorithm compute the slope between the two dose distributions [14].