5. Application to 3D microwave radar images: Measuring changes within a reconstruction (Case 1) - djkurran/Automated-framework-for-evaluating-microwave-and-multi-modality-breast-images GitHub Wiki
The first of two examples is presented to demonstrate the application of the workflow to three-dimension (3D) backscattered energy images; the previous examples applied the segmentation and analysis workflow to 2D scenarios. Moreover, for the examples presented in sections 2-4, a dual-modality breast imaging approach was taken by combining ultrasound and microwaves. The workflow was used as a tool to test a more general hypothesis that incorporating prior information extracted from ultrasound imagery into a reconstruction algorithm improves the quality of the images reconstructed relative to present state-of-art variants. In section 2 the refinements were related to the use of prior information to construct an inhomogeneous numerical background that is used by the reconstruction algorithm. In sections 3 and 4, the change was related to a refinement of the path ray model used for the time-delay calculations.
For the examples presented in sections 5 and 6, the reconstruction algorithm is fixed over the numerical experiments that are conducted. The objective is to demonstrate how the workflow 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. Accordingly, the quantitative results may be used to support or challenge an hypothesis held by the researcher.
The image processing and analysis workflow presented in sections 2 – 4 has been adapted to operate on volumes and 3D representations of masks (or objects). Furthermore, these 3D microwave radar examples differ from the 2D examples in that the aim is to compare and quantify differences between two reconstructions. For this scenario, the parameters that are initially implemented by the reconstruction operator are set to base case values. A reference image is formed when the operator is applied to microwave reflection data.
The parameter values used by the operator are subsequently modified, and the operator is applied to the same microwave reflection data. This leads to a reconstruction that is referred to as a test image. Application of the segmentation and analysis workflow to the reference and test images is used to measure changes within the test image that arise due to a variation of a reconstruction operator parameter.
Figure 1. Image reconstruction workflow [8, 23, 24].
A general overview of the 3D reconstruction workflow and the reconstruction operator that is applied to the backscattered fields is presented in section 5.1. The overview provides key background information and identifies important parameters used in the image reconstruction process. The numerical experiment used to demonstrate how the segmentation and analysis workflow may be used to analyse results is described in section 5.2. The segmentation methodology that is applied to the reference and test images is presented in section 5.3. The aim of the methodology is to delineate regions associated with reflections from significant scatterers. A general backscatter energy analysis of these dominant responses present within the test and reference images is carried out and the results are shown in section 5.4. The metrics described in section 1.4.b(1) have been adapted for the application to volumes and 3D objects, and so are applied to the test and reference 3D masks in order to evaluate the geometric properties of the dominant scattering regions within the reconstructions. The geometric analysis results are presented in section 5.5. Measurable changes in the response of a dominant scatterer, smearing of a response, changes in the intensity and focus of the response, and changes to the extent of a response are examined in section 5.6. Finally, in order to provide an understanding beyond the dominant responses, the structural similarity index measure (SSIM) and normalized root mean square error (NRMSE) metric are applied over 2D slices along each of the coordinate axes. The results are presented in section 5.7.
Figure 2. Scattered density breast model constructed from an MRI scan of a patient. Gland consists of three different types of tissues indicated by different colors. A
Table I. Scattered density breast model Debye parameters [21]
A general overview of the workflow implemented to reconstruct 3D images from microwave radar data is provided. The workflow is summarized in Figure 1.
The first step of the workflow is to generate a set of data. For these examples, data are generated by simulating a radar system operating in the monostatic mode. In this mode, a single balanced antipodal Vivaldi antenna with a director included in the aperture (BAVA-D) [20] illuminates a breast model with an ultra-wideband pulse. The 3D numerical model of a scattered density breast is described in [21] and is constructed from an MRI scan acquired from a patient study reported in [22]. The model is covered with a
The model and antenna are immersed in Canola oil and are constructed with SEMCAD X (SPEAG, Zurich). The setup is simulated with a finite difference time-domain (FDTD) solver which computes the electromagnetic fields. The single antenna is moved to
The backscattered field computed at each antenna position is extracted by subtracting the incident field from the total field. The total and incident fields correspond to the fields computed in the presence and absence of the numerical breast model, respectively. Once the data are generated, an adaptive filtering technique is applied to each backscattered field to remove the dominant response corresponding to the skin reflection [23]. This results in the set of preprocessed time-domain signals, { s1(t), s2(t) , ..., sM(t) }.
To form an image, the set of preprocessed signals are synthetically focused to a specific location, or focal point, in the imaging domain that represents the breast interior. A backscattered signal that originates from a focal point propagates along a path ray to arrive at the synthetic array with different times-of-arrival at each antenna location. The delay-and-sum image reconstruction operator compensates for these delays by applying a reverse delay to each preprocessed signal. This has the effect of time shifting the signals to align the reflections from the focal point. Accordingly, if the time delays are accurately computed, the data are focused so that the signals from each antenna add constructively at the focal point.
The time delays are computed with two preprocessing steps. First, the breast surface is reconstructed. For the numerical examples, a tessellated representation of the breast surface is extracted from the numerical model. The surface is used to construct the imaging domain that is comprised of voxels bound by the surface [24].
For the second step, the path ray that extend from a focal point at
Next, real valued permittivities are assigned to each segment of the path ray; the material within each partition is assumed to be homogeneous, lossless, and constant across the band of frequencies that comprise the ultra-wideband signal (i.e., non-dispersive). This information is used to evaluate the velocity of the signal over the path length. Finally, the round trip travel time of the signal over the path is computed using the segment velocities and lengths. For each focal point, a set of
During focusing, the reconstruction operator applies the corresponding time-delay,
The focal point iteratively moves to each reconstruction location within the breast, so the process is repeated over all voxels within the imaging domain. The result is a 3D backscattered energy image of the breast interior. Importantly, all time-shifted responses that arise from significant scatterers, such as regions of malignant tissue, are coherently summed; reflections from healthy tissues and any artifacts from preprocessing are assumed to add incoherently. Consequently, regions within the image where the backscatter energy intensity is elevated may be interpreted as responses corresponding to dominant scattering regions. These responses may, in turn, be used to detect the presence of a malignancy.
Figure 3. Reference image. A permittivity of
Figure 4. Test image. A permittivity of
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 (the intensity of a response, for example). A key step in microwave radar image reconstruction is the time delay calculation. The calculation is dependent on a number of variables that may impact the accuracy of the calculation. 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.
As described in section 5.1, an image is formed by synthetically focusing signals to focal points within the imaging domain. A signal that originates from a focal point, travels along a path ray and arrives at the synthetic array with different times-of-arrival at each antenna location. A path ray is modeled with three constituent segments: immersion medium, skin, and breast interior. The length of each segment varies and is dependent on the focal point location and direction (or angle) of the ray [24].
Real valued permittivities are assigned to each segment to evaluate the velocity of the signal over the path ray. The round trip travel time of the signal over the path is computed using the segment velocities and lengths. These travel times, or time delays, are computed by assigning
However, to form the reference image, a value of
The relative permittivity values are used in the preprocessing step described in section 5.1 to compute the time-delays for all path rays. 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
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. Namely, the value assigned to all breast interior segments to form the test image is increased by
Figure 5a. Segmentation and analysis workflow.
Figure 5b. Overview of automated workflow to segment and evaluate 3D microwave radar breast images
The reference and test images shown in Figures 3 and 4, respectively, are input to the workflow represented in Figure 5. As described in section 5.2, the images are comprised of voxels that are interior to the breast surface. The entire image corresponds to the region-of-interest, so no pre-processing of the image is carried out prior to segmentation.
The unsupervised machine learning based segmentation algorithm is applied to the reference and test 3D images (not individual slices within the 3D images) to partition the volumes based on the intensity of the reconstructed backscatter energy. The algorithm applies the
When the segmentation algorithm is terminated, the reference and test volumes are partitioned into
Figure 6. Reference image segmentation results. Reference image (left column) is input to segmentation workflow and partitioned into
Figure 7. Test image segmentation results. Test image (left column) is input to segmentation workflow and partitioned into
The segmentation algorithm partitions the 3D image based on the value of the intensity of the reconstructed backscatter energy. The cluster-to-mask mapping function is simplified to suit the functional interpretation of the backscatter energy image. Given the region-of-interest has been partitioned into
The construction of masks leads to a set of test and reference objects, where object refers to a region of connected voxels within the mask. 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 has 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 8. 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 9. 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.
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,
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 8.
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,
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 9.
Figure 10. Detailed geometric analysis to evaluate the similarity between
As described in section 5.2, the reconstruction operator used to form the reference and test images is applied to the same backscattered fields calculated from the equivalent forward model. The backscatter energy intensities are reconstructed onto 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.
A comprehensive geometric property analysis is performed to assess the overlap between a test object under investigation and the connected reference object using the metrics described in section 1 that include metrics 1-8 (i.e., Dice Coefficient, Jaccard, Ratio-of-detection, Artefact-rejection). The similarity in shape is assessed with metric 16, the average Hausdorff distance. If the test object is not connected with any of the other reference objects, then the analysis is not conducted.
For the set of segmented test objects, test object
Recall from section 5.2 that the reference image is set to an image reconstructed by an operator with a baseline set of parameter values. The parameter values assigned to the reconstruction operator that was applied to the backscattered fields to reconstruct the test image differ from the baseline values by
The results are not included in the submitted manuscript [11], as the overlap between the test and reference object is not significant. As pointed out in section 1.4, the metric values are most meaningful when the ground truth is used as the reference image. That is, when the responses reconstructed within the test image are compared with malignant tissue regions within the forward model. For scenarios where differences between reconstructions are being evaluated, the geometric analysis results may be more relevant for scenarios where a controlled variable is perturbed by a small amount.
Figure 11. Comparative backscattered energy analysis of
The backscatter energy intensities within the test object and the nearest reference object are compared. 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
Starting with the first test object,
Likewise, the maximum intensity,
The test image is formed with a reconstruction algorithm that uses permittivity values
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.
Conversely, 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
Differences between the test and reference response are evaluated with metrics
Next, the volume ratio given by (29) measures the ratio between volumes
The change in FWHM,
Finally, the change in SMR relative to the reference,
Collectively, the comparative analysis infers that the increase in independent variable value leads to a decrease in accuracy of the calculated time-delays which has a negative impact on image quality. Although the response is sharper and more focused, the reconstructed response has shifted, is weaker, and there is an increase in clutter within the imaging domain. The results are summarized in Figure 11.
The comparative analysis is repeated for the remaining test objects. The results for
Figure 12. Comparative backscattered energy analysis of
Figure 13. Comparative backscattered energy analysis of
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
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 a researcher 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 along each point 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 14, 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 <
The forward model tumor is located at
To complement the analysis using the NRMSE metric, the Structural Similarity Index (SSIM) presented in [10] 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
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 14 and can be compared with the corresponding NRMSE plot shown in the bottom row. The SSIM plots do not appear to furnish meaningful information related to changes within the test image compared to the NMRSE plots. As presented in [10], the motivation for the development of the metric is to quantify subtle changes to image quality 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. 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 14. 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.