6. Application to 3D microwave radar images: Measuring changes within reconstructions (Case 2) - djkurran/Automated-framework-for-evaluating-microwave-and-multi-modality-breast-images GitHub Wiki


The second example is presented to demonstrate the application of the segmentation and analysis framework to three-dimension (3D) backscatter energy images. Similar to the first example, the objective is to demonstrate how the framework may be used assist a researcher who is conducting an experiment to investigate how a change to an independent variable used in the reconstruction process impacts a dependent variable associated with a response that is reconstructed within the image.

The breast model, numerical data and reference image are the same as the first example. Parameter values used by the operator are revised compared to the base line case, and the operator is applied to the same microwave reflection data to form a test image. Application of the segmentation and analysis framework to the reference and test images is used to measure changes within the test image that arise due to variations of the parameter implemented by the reconstruction operator.

Background information and details related to the numerical experiment used to demonstrate how the segmentation and analysis workflow may be used to analyse results is described in section 6.1. The segmentation methodology is applied to the reference and test images and the results are presented in section 6.2. A general backscatter energy analysis of the dominant responses detected within the test and reference images is carried out and the results are presented in section 6.3. The geometric analysis of the dominant scattering regions within the reconstructions is presented in section 6.4. A detailed analysis is carried out to compare each response reconstructed in the test image to the corresponding response reconstructed in the reference image. The results are presented in section 6.5. The normalized root mean square error (NRMSE) metric and the structural similarity index measure (SSIM) are applied over 2D slices along each of the coordinate axes. The results are presented in section 6.6.

Figure 1. Reference image. A permittivity of $\epsilon_{r,int} = 4.47$ is assigned to all interior segments of the path rays when computing the corresponding time-delays. The reconstruction operator applies these time-delays to the preprocessed backscatter signals when forming an image. The reconstruction is sliced along each of the coordinate axes at the location of maximum intensity.

Figure 2. Test image. A permittivity of $\epsilon_{r,int} = 16.0$ is assigned to all interior segments of the path rays when computing the corresponding time-delays. These values deviate from the values assigned to the interior path ray segments used to compute the time-delays for the reference image by $258$ %. The reconstruction operator applies these time-delays to the backscatter signals when forming an image. The image is sliced along each of the coordinate axes at the location of maximum intensity.

6.1 Description of Case 2


This case continues with the first example presented in section 5. The objective is to demonstrate how the segmentation and analysis framework may assist a researcher who is investigating how a change to an independent variable such as a reconstruction operator parameter impacts a dependent variable. Similar to the first example, an experiment is conducted in which an independent variable used in the time delay calculation is changed relative to a reference value. The segmentation and analysis workflow is applied to the test and reference images to evaluate the impact the change has on a dependent variable, such as the intensity of a response.

The time delays are computed by assigning $\epsilon_{r,ext}=2.5$, and $\epsilon_{r,skin}=36.4$, to the immersion medium and skin components of the ray path, respectively. The same immersion medium and skin segment value assignment is implemented for both the reference and test images.

Just as in the first case, to form the reference image, a value of $\epsilon_{r, int}=4.47$, is assigned to the breast interior segments. Conversely, to form the test image, a value of $\epsilon_{r, int}=16.0$ is assigned to the interior segments. This means that the value assigned to all breast interior segments for the test image case is $258$ % higher compared to the value used for reference case.

The relative permittivity values are used in the reconstruction preprocessing step described in section 5.1 to compute the time-delays for all path rays. Recall that during focusing, the reconstruction operator applies these time-delays to the preprocessed signal to time-shift each signal by the amount specified by the time-delay. The operator then combines the time-shifted signals over all $M$ antenna positions of the synthetic array and the resulting value is assigned to the voxel. The value assigned to the voxel represents the intensity of the backscatter energy at the focal point within the imaging domain. The resulting backscatter energy images assigned to the reference and test images are shown in Figures 1 and 2, respectively.

In the context of a controlled experiment, the independent variable is the relative permittivity assigned to the interior segments of the path rays that are used to compute the time delays. The segmentation algorithm delineates dominant scattering regions that are possibly associated with malignant tissue. The analysis framework is applied to the segmented regions to measure the affect that the change of the independent variable has on a dominant scattering region. Measurable changes of a dependent variable include shifts in the response of the dominant scatterer, smearing of a response, changes in the intensity and focus of the response, and changes in the extent of a response.

Figure 3. Reference image segmentation results. Reference image (left column) is input to segmentation workflow and partitioned into $12$ clusters (right column). The clusters are mapped to a tissue type image (center column). Figure 4. Test image segmentation results. Test image (left column) is input to segmentation workflow and partitioned into $13$ clusters (right column). The clusters are mapped to a tissue type image (center column).

6.2 Segmentation of reference and test images


The reference and test images shown in Figures 1 and 2, respectively, are input to the segmentation and image analysis workflow that is described in section 5.2. The images are comprised of voxels that are interior to the breast surface, so the segmentation algorithm is applied directly to the imaging domain.

Similar to the 3D microwave radar example presented in section 5, the unsupervised machine learning based segmentation algorithm is applied to the reference and test 3D images to partition the volumes based on the intensity of the reconstructed backscatter energy. The algorithm applies the $k$-means clustering algorithm to the image in an iterative fashion. After each iteration, a statistical test evaluates the hypothesis that the cumulative distribution function of the voxel values within clusters at the previous and current iterations originate from different distributions. If the hypothesis does not fail, then the number of clusters used by the algorithm is incremented, otherwise convergence is implied and the volume is not partitioned further. The iterative procedure is summarized by the flow diagram presented in Figure 5 of 2.1.3.

When the segmentation algorithm is terminated, the reference and test volumes are partitioned into $12$ and $13$ clusters, respectively. The segmentation results where the 3D image is sliced at the location of maximum intensity sliced along each of the coordinate axes are shown for the reference and test volumes in Figures 3 and 4, respectively. The segmentation algorithm partitions the 3D image based on the intensity of the reconstructed backscatter energy. The cluster-to-mask mapping function maps the highest valued clusters to voxels within a binary image comprised of significant scattering masks. The remaining clusters are mapped to voxels within another binary image related to the region for which there is not significant scattering.

The construction of masks leads to a set of test and reference objects. For some scenarios, the presence of a large number of significant scattering regions related to glandular/fatty tissue interfaces may occur. This complicates the application of the 3D metrics, as some of these regions may be considered as image noise (e.g. small relative to the breast interior). Accordingly, the post-processing technique developed for the microwave tomography images been adapted to filter out these small isolated regions or artefacts. The volume enclosed by each closed surface is evaluated, and the ratio of each of the volumes to the largest volume is computed. Ratios below a user-defined extraction threshold are removed from the analysis.

In order to enhance the flexibility to identify artefacts and regions of interest, all significant scattering regions are automatically labeled as part of the mask construction process. This allows the user to specify which of the labeled scattering regions to remove from further analysis (small objects near the skin surface or chest wall, for example).

An analysis framework comprised of metrics described in section 1 is applied to the final set of test and reference objects. The aim of the analysis is to measure the affect that the change of the reconstruction operator parameter has on the response of the dominant scatterer that has been segmented from the reference and test images. The analysis is described next.

Figure 5. Dominant scattering regions segmented from test image that result in a set of test objects. The objects are displayed within the tessellated breast surface that bounds the imaging domain. The results of a general analysis of the backscatter energy associated with each of the detected dominant scattering regions is shown in the table.

Figure 6. Dominant scattering regions segmented from reference image that result in a set of reference objects. The objects are displayed within the tessellated breast surface that bounds the imaging domain. The results of a general analysis of the backscatter energy associated with each of the detected dominant scattering regions is shown in the table.

6.3 General backscatter energy analysis


A general analysis is carried out by applying backscatter energy analysis metrics, described in section 1.4.b(1), to the test and reference objects. Each test object is applied to the test image to extract voxels with values corresponding to the backscattered energy intensity within a region associated with a dominant scatterer. The maximum intensity, $I_{max,test(i)}$, over all voxels within the test object is determined with (23). The coordinates of the location of this value, $I_{max,test(i)}(\vec{r})$, are determined and used to calculate the FWHM and the SMR of the response using (21) and (22), respectively. Finally, the volume of each test object is evaluated with (24). The metrics that are applied to each object are described in detail in section 1.4.b(1).

The maximum backscattered energy intensity of each test object is scaled to the maximum value over the imaging domain of the test image to assist with the interpretation of the results. The value, along with the other values calculated for the general analysis are presented in a table. The objects are displayed within the imaging domain that is bound by the tessellated breast surface. The objects and tables are shown in Figure 5. The analysis is repeated for each of the reference objects. That is, each reference object is applied to the reference image to extract voxels with values corresponding to the backscatter energy intensity within a region associated with a dominant scatterer. The maximum intensity, $I_{max,ref(j)}$ over all voxels within the reference object is determined with (25). The coordinates of the location of this value, $I_{max,ref(j)}(\vec{r})$, are determined and used to calculate the FWHM and the SMR associated with the response. The volume of each reference object is calculated with (26).

To assist with a qualitative comparison of results, a similar table and display are provided for the reference objects. Note that the maximum backscattered energy intensity that is displayed for reference object is scaled to the maximum value over the imaging domain of the reference image. The objects and tables are shown in Figure 6.

6.4 Geometric property analysis


As described in section 5.2 for the first example, the backscatter energy intensities for the test and reference images are reconstructed on voxels within an imaging domain with the same coordinates. Therefore, as a preprocessing step, the reference objects are transformed to a space shared by the test object.

The test object under investigation is then examined to determine if it is connected to any of the reference objects. The objects are connected if they share any of the same voxels. If this is the case, then the test and connected reference objects are evaluated within a region of interest defined as the box that bounds the union between the reference and test objects.

None of the test objects are connected to a reference object, so a detailed geometric analysis is not conducted on these objects; no results are reported.


Figure 7. Comparative backscattered energy analysis of $test(1)$. It is compared with the nearest reference object, $ref(1)$. Deviation between location of responses is $38$ $mm$.

Figure 8. Comparative backscattered energy analysis of $test(2)$. It is compared with the nearest reference object, $ref(1)$. Deviation between location of responses is $30$ $mm$.

Figure 9. Comparative backscattered energy analysis of $test(3)$. It is compared with the nearest reference object, $ref(3)$. Deviation between location of responses is $24$ $mm$.

6.5 Comparative backscatter energy analysis


The backscatter energy intensity associated with a test object is compared with the nearest reference object. Similar to the preprocessing step described in section 5.5, the reference objects are transformed to a space shared by the test object. For each test object, the nearest reference object is identified (i.e., shortest distance between centroids). For example, for the $i^{th}$ test object $test(i)$ the nearest reference object is the $j^{th}$ reference object $ref(j)$.

Starting with the first test object, $test(1)$, the nearest reference object is $ref(1)$. The maximum intensity, $I_{max,test(1)}$, over all voxels within the test object is determined with (23). The coordinates of the location of maximum value within the test object are determined and are used to calculate the FWHM and the SMR of the response using (21) and (22), respectively.

Likewise, the maximum intensity, $I_{max,ref(2)}$ over all voxels within the reference object $ref(1)$ is determined with (25). The coordinates of $I_{max,ref(2)}$ are used to evaluate the shift in the test response relative to the reference image with (27). Differences between the test and reference response are evaluated with metrics 28−31. The comparative analysis results for $test(1)$ are presented in Figure 7.

The comparative backscatter energy analysis is repeated for test objects $2$ - $3$. The results for $test(2)$ and $test(3)$ are shown in Figures 8 and 9, respectively. For the first three test objects examined, the results presented in Figures 7-8 demonstrate that the nearest reference object used for the comparative analysis is a significant distance from the test object under investigation. This is quantified by the shift-in-response values measured for $test(1)$, $test(2)$, and $test(3)$, which are $38$ $mm$, $30$ $mm$, and $24$ $mm$, respectively. For these test objects, the nearest reference object is not in close proximity. Furthermore, none of these test objects corresponds to the dominant response within the test image. It is unlikely that any of the first three test responses are compared with a corresponding reference response.

The analysis of the remaining test object, $test(4)$, is of more interest compared to the first three test objects as the dominant response within the test image corresponds to $test(4)$. Moreover, the dominant response within the reference image corresponds to the nearest reference object, $ref(2)$. The maximum intensity, $I_{max,test(4)}$ over all voxels within the test object is determined with (23). The coordinates of the location of maximum value within the test object, $I_{max,test(4)}(\vec{r})$, are < $-18, -10, -12$ > $mm$ and are used to calculate the FWHM and the SMR of the response using (21) and (22), respectively.

Likewise, the maximum intensity, $I_{max,ref(2)}$ over all voxels within the reference object $ref(2)$ is determined with (25). The coordinates of $I_{max,ref(2)}(\vec{r})$ are $<-20, 2, -28>mm$ and used to evaluate the shift in the test response relative to the reference image with (27).

The test image was reconstructed using time-delays calculated by setting the relative permittivity of the interior path ray segments, $\epsilon_{r,int}$, to a value that deviated from the mean value by $258$ %. The shift-in-response that is calculated indicates that the test response has shifted by $20 mm$, relative to the reference response. As already noted, the centroid of the tumor represented in the forward model tumor is $<-20, 0, -28>mm$. Therefore, the maximum intensity, $I_{max,test(4)}$ is displaced from the centroid of the ground truth tumor implemented by the forward model by $19$ $mm$.

For comparison, the reference image was reconstructed using time-delays calculated by setting the interior path ray segments to the mean relative permittivity calculated over all voxels of the breast model. The dominant response within the image deviated from the ground truth tumor centroid by $2$ $mm$. The test image used in the first example was reconstructed using time-delays calculated by setting the relative permittivity of the interior path ray segments, $\epsilon_{r,int}$, to a value that deviated from the mean value by $57$ %. The dominant response within the image deviated from the ground truth tumor centroid by $7$ $mm$, and the reference response by $9$ $mm$. The shift-in-response from the reference image and ground truth are displayed in Table 1.

When accurate time-delay values are used, the reconstruction operator is applied to signals that are correctly time-shifted. The time-shifted signals are combined to form an image, and corresponding peaks in each signal interfere constructively where scatterers are located; the signals combine destructively in other regions of the image.

The test images reconstructed for case 1 and 2 used relative permittivity values that deviated by $57$ % and $258$ % from the value used to reconstruct the reference case. As part of the reconstruction workflow (see Figure 1 as part of section 5.1), the algorithm assigns the values to interior segments of each path ray connecting an array element to a focal point. The interior segment path values were used to calculate the time-delays implemented by the reconstruction operator to form the test image. Time-delay values that deviate from actual values perturb the amount of time-shift applied to each preprocessed signal.

When the reconstruction operator is applied to signals with perturbed time-shifts that arise from inaccurate time-delays, signal peaks no longer center on the scatterers. Consequently, constructive interference of signals occurs in locations that are offset from the actual scatterer locations. This offset in constructive interference manifests as a shift in the test response relative to the centroid of the actual tumor. An increase in the shift in the response is observed (refer to Table 1) as the time-delay calculation deteriorates.

Table 1. Comparing backscattered energy between dominant scatterer within reference, test case 1, and test case 2 images. $\epsilon_{r,ext}$ and $\Delta \epsilon_{r,ext}$ % represent the value set to the relative permittivity of the interior segment of the ray path and the change in value from the reference case, respectively.

Differences between the test and reference response are further evaluated with metrics 28−31. First, the ratio between the maximum intensity within $test(4)$ and the maximum intensity within $ref(2)$ and is evaluated with (28). The resulting intensity ratio is $0.48$ which implies that the response reconstructed in the test image is diminished, relative to the response in the reference image. For comparison, the intensity ratio of the corresponding maximum response within the test image of the first example was $0.64$, and is shown Table 1.

Perturbations of time-shifted signals due to the use of inaccurate time-delay leads to a misalignments of signal peaks that arise from tissue interfaces. When the signals are combined, this misalignment manifests as a diminished and weaker response relative to a response reconstructed from accurate time-delays. The diminished intensity of the dominant response reconstructed in the test images relative to the reference image supports the observation that the accuracy of the time-delay calculation deteriorates as the value, $\epsilon_{r,int}$, assigned to the interior segments of the path ray deviates from the value implemented to reconstruct the reference image.

Next, the volume ratio given by (29) measures the ratio between volumes $test(4)$ and the nearest reference object $ref(2)$. The volume ratio is intended to complement the FWHM metric. The volume ratio is $0.42$, which implies that the test response is more compact, focused, and sharper relative to the reference image. A more compact response relative to the reference response is observed qualitatively by comparing $test(4)$ and $ref(2)$ shown in Figure 10.

The change in FWHM, ΔFWHM, calculates the difference in FWHM between the $test(1)$ and nearest reference object $ref(2)$, and is given by (30). A negative distance of $−4$ $mm$ is measured which suggests the extent of the test response is reduced and more compact relative to the extent of the reference response. The result supports the reduced volume ratio value.

Finally, the change in SMR relative to the reference, $\Delta SMR$ %, is calculated with (31) to be $−7.7$ %. The negative $\Delta SMR$ % is attributed to the receding response intensity. However, the negative impact of the diminished intensity is countered to some extent by reduced background clutter and artifacts external to the response. This is counter to what was implied for Case 1 where the decrease in SMR is attributed to a decrease in response intensity and an increase in image clutter.

Collectively, the comparative analysis of $test(4)$ with the corresponding reference infers that the increase in independent variable value leads to a decrease in accuracy of the calculated time-delays which has a significant negative impact on the reconstruction of the dominant scatterer. Similar to Case 1 described in section 5, although the response is sharper and more focused, the reconstructed response has notably shifted, and is weaker. However, an improvement in image quality is implied by a decrease in clutter and artifacts relative to both the reference and Case 1. This is supported qualitatively when the cluster images shown in Figure 3 and Figure 7 of section 5 are examined. Namely, significant scattering (or artifacts) are observed near the chest wall for Case 1 and the reference, but is not as significant for the Case 2. The results are summarized in Figure 10.

Figure 10. Comparative backscattered energy analysis of $test(4)$ which corresponds to the dominant scatterer within the test image. It is compared with the nearest reference object, $ref(2)$ which corresponds to the dominant scatterer within the reference image. The dominant response within the test image has been shifted by $20$ $mm$ relative to the reference image.


6.6 Comparative image analysis of 2D slices over each coordinate axis

Two-dimensional slices of the test and reference imaging domains are compared along each of the coordinate axes with the Normalized Root Mean Square Error (NRMSE). The metric is given by (17) except $ref_{mask}$ and $test_{mask}$ are replaced with two-dimensional slices of the reference and test imaging domains, respectively, at a point along a coordinate axis.

As described in section 1.4b(4), each slice is preprocessed by removing pixels not contained in the imaging domain. This step is taken, as inclusion of these background pixels may reduce the accuracy of the metric. The metric is then applied to vectorized versions of the preprocessed slices to measure the normalized distance between the test and reference vectors. Hence, it heuristically informs us how close the test image is to the reference image. A value of zero indicates a perfect match between images; an increasing value indicates that the difference between test and reference images is increasing. For the 3D examples, the NRMSE implies how much the test image has changed relative to the reference image due to the perturbation of the reconstruction operator parameters.

The metric is applied to slices at equally spaced points along a coordinate axis resulting in a set of values that are then plotted. The process is repeated for each coordinate axis. The plots are shown along the bottom row of Figure 11, and may be used as a tool to evaluate how the test and reference images deviate from each other along each of the coordinate axes.

To assist with the interpretation of the plots, note that the coordinates of the location of maximum value within the test and reference images is < $-18,−10,−12$ > $mm$ and < $−20,2,−28$ > $mm$, respectively. The centroid of the forward model tumor is < $−20,0,−28$ > $mm$.

The forward model tumor is at $y=0$ $mm$ and $z=-28$ $mm$, and the axial and coronal plots show moderate NRMSE values in this region. This suggests some difference in these images such as a shift in the location of the responses, as expected in Figures 1, 2, and 10. In the sagittal plane, the NRMSE is lowest at $x=−18$ $mm$ which is in close proximity to the physical tumor location and the maximum intensity of the reference image in this plane.

To complement the analysis using the NRMSE metric, the Structural Similarity Index (SSIM) presented in [6] is applied to two-dimensional slices of the test and reference imaging domains along each of the coordinate axes. It measures perceptual differences between two similar images. The computed values range from 0 to 1, where 1 indicates that there are no perceptual difference between the test and the reference images.

Similar to the NRMSE metric, each slice is preprocessed by removing pixels not contained in the imaging domain. The SSIM Structural Similarity Index operator given by (32) is then applied to the preprocessed slices along each of the coordinate axes. A more detailed description is provided in section 1.4b(4).

The plots are presented in top row of Figure 11 and can be compared with the corresponding NRMSE plot shown in the bottom row. As presented in [6], the motivation for the development of the metric was to quantify image quality degradation caused by processing such as data compression or by losses in data transmission.

For the 3D reconstructions presented in sections 5 and 6, the changes in the dependent variable are large relative to the base line case ( $57$ %, $258$ %, respectively). These perturbations lead to significant spatial shifts, variations of the extent, and distortions of shape and intensity of the reconstructed responses. The dependent variable changes also manifest as image artefacts. These reconstructed response changes differ from the small subtle changes that the metric was intended to quantify. Hence, other metrics such as the NRMSE appear to be better suited for evaluating image changes compared to the SSIM. Nevertheless, it has been as a tool for assessing image differences by microwave imaging researchers (see [7, 26], for examples).

Figure 11. SSIM applied to series of 2D slices (bottom row) along axial (left), coronal (center), and sagittal direction (right) for comparison between test and reference images. Application of NRSME to series of 2D slices shown on top row.