Iterative peak combination: a robust technique for identifying relevant features in medical image histograms

Histogram-based methods can be used to analyse and transform medical images. Histogram specification is one such method which has been widely used to transform the histograms of cone beam CT (CBCT) images to match those of corresponding CT images. However, when the derived transformation is applied to the CBCT image pixels, significant artefacts can be produced. We propose the iterative peak combination algorithm, a novel and robust method for automatically identifying relevant features in medical image histograms. The procedure is conceptually simple and can be applied equally well to both CT and CBCT image histograms. We also demonstrate how iterative peak combination can be used to transform CBCT images in such as way as to improve the Hounsfield Unit (HU) calibration of CBCT image pixel values, without introducing additional artefacts. We analyse 36 pelvis CBCT images and show that the average difference in fat tissue pixel values between CT images and CBCT images processed using the iterative peak combination algorithm is 23.7 HU. Compared to 136.7 HU in unprocessed CBCT images and 50.9 in CBCT images processed using histogram specification.


Introduction
The question we address in this work is one of transforming a cone beam CT (CBCT) (Jaffray et al 2002, Létourneau et al 2005 image histogram to more closely resemble that of a corresponding planning CT image histogram. This is a commonlyoccurring problem, given the well-studied artefacts (Siewerdsen and Jaffray 1999, 2001, Poludniowski et al 2011 in CBCT images which often result in image pixel values that do not correspond to Hounsfield Units (HU) (Fotina et al 2012).
Histogram-based methods are widely-used tools for analysing and processing medical images and have been suggested as a method to improve HU calibration in CBCT images. A typical method for addressing this problem is the use of histogram specification, also known as histogram matching, algorithms which aim to transform one image so that its histogram is the same as another. Cumulative frequency distributions are calculated from image histograms and used to determine how to adjust the shape of the source image (CBCT) histogram to match the reference image (CT) histogram. Amit and Purdie (2015), for example, use histogram specification to enable CBCT images to be used in a procedure which automatically generates radiotherapy treatment plans for breast cancer patients; Park et al (2017) use histogram specification to improve the performance of CT-CBCT deformable image registration; Zhang et al (2017) use histogram matching to allow CBCT images to be used for radiotherapy dose calculations; and Arai et al (2017) use histogram specification to enable CBCT images to be used for proton therapy dose calculations.
The main deficiencies of the histogram specification algorithm when used in this context are that CT image histograms contain far more detail than CBCT histograms, and that there may be anatomical changes in the patient between the acquisition of the CT and CBCT images. When histogram specification is used to transform a CBCT image histogram to match a CT histogram, detail must be added to the CBCT histogram. However there is no correspondence between the voxels entering a bin of a histogram and their Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. spatial location in the image. Therefore the detail that is added, for example increased contrast between fat tissue and muscle tissue, is unlikely to correspond to the location in the CBCT image where the detail is required. Anatomical changes in the patient between CT and CBCT imaging (e.g. weight loss), or different fields of view, will result in different numbers of voxels representing each tissue or material being present in each of the images. Therefore the true CBCT histogram is not expected to exactly match that of the CT image. Forcing the histogram of one image to match the other will risk artificially changing one material to another in some voxels.
We approach the problem from the opposite direction; using the CBCT image as the target and distilling the information in the CT histogram until it resembles the CBCT. It then becomes more straightforward to compare the CT and CBCT histograms, and a more simple transformation can be defined which maps the CBCT image histogram to the CT. This approach has previously been taken by Marchant et al (2008), where the CT image histogram was simply smoothed until sufficient detail had been removed, and by Thing et al (2016), who used histograms with wide bins to reduce the amount of detail and enable the identification of peaks in the histograms corresponding to simply soft tissue and air. However, it is likely that the amount of smoothing or number of histogram bins used would need to be adjusted based on the particular patient, the anatomical site of the image or the acquisition parameters of the scanner.
In view of the recent interest in image correction methods based on histogram features, particularly in the field of CBCT imaging, it is important to find improved methods to extract features and map one image histogram to another. Therefore we propose an automatic method for identifying relevant features present in histograms of CBCT and corresponding planning CT images. The iterative peak combination algorithm has a single adjustable parameter and is conceptually very simple. We then use the results of the iterative peak combination algorithm to determine a linear transformation which scales pelvis CBCT image histograms to match the corresponding CT histograms. The same transformation is then applied to the CBCT image pixel values. An automated analysis of the raw CBCT image pixel values is performed and compared to results obtained after using iterative peak combination and histogram specification to process the CBCT images.

Image histograms
The image histograms created in this work are comprised of bins, each with a bin width, and defined by two bin edges e i low and e i high . The height of each bin corresponds to the frequency, f i , of image pixels falling into that bin. For a given image, the image histogram is created with a chosen bin width. The number of histogram bins is determined by finding the range of pixel values present in the image and dividing by the bin width. Every pixel in the image is then looped over to calculate which bin the pixel value falls into and the frequency of that bin is incremented by one.
The first bin of the histogram typically contains entries corresponding to pixels which fall outside the field-of-view of the scanner that acquired the image, and which all have the same value. The frequency of the first bin is therefore automatically set to zero before continuing.

Histogram peak identification
A peak, p, in the histogram is defined as being located at the point at the centre of a bin that has frequency larger than both neighbouring bins, i.e. This definition is clearly susceptible to finding spurious peaks due to noise in the histogram. This problem can be mitigated by using a larger bin width when creating the image histogram. However, the iterative peak combination technique described below is robust against peaks produced by noise. The use of wider histogram bins, or additional histogram smoothing, is therefore not required.
Once all peaks in the histogram have been identified, the peaks are iteratively combined using the following procedure.

Iterative peak combination
Iterative peak combination is inspired by the problem of jet-finding in high-energy particle physics (Salam 2010). The theory underlying jet-finding is very well-motivated for the applications in particle physics, and is not relevant to the discussion here. However the general idea of the technique-transforming a large amount of information into a small number of salient features-is simple and easily adapted to the problem of locating relevant peaks in histograms of medical images.
The peak combination algorithm proceeds as follows: (i) Calculate the distance d i j , between all peaks p i and p j (see equation (2)).
(iii) Combine the two closest peaks p i min and p j min to form peak p ij (see equation (3)).
(iv) Remove peaks p i min and p j min from the list of peaks.
(v) Add combined peak p ij to list of peaks.
(vi) Repeat until number of remaining peaks equals n target .
The distance d i j , is defined as where p i x , and p i y , are the location (in HU) and height, respectively, of peak i.
The location of the combined peak, p ij x , , is constructed by taking the average of the locations of peaks i and j, weighted by the height of the peaks. The height p ij y , is calculated by linearly interpolating between peaks i and j at the location p ij x , , The histogram-scaled CBCT image is then produced by applying the transformation T to each pixel in the image.

Histogram specification
Histogram specification, also known as histogram matching, is a widely-used technique for transforming the histogram of one image to match that of a specified reference image (Gonzales and Woods 2008). The cumulative frequency distributions of the CT and CBCT histograms are used to calculate a transformation that transforms a given CBCT histogram bin frequency to match the corresponding bin of the CT histogram. We refer to the resulting image as the histogram-specified image.
In this study we use the implementation of histogram specification provided in the Insight Toolk-it (Yoo et al 2002).

Histogram-transformed image comparison
CT and CBCT images of patients undergoing radiotherapy of the pelvis area are used to demonstrate and test the histogram transformation techniques. The pixel values in the histogram-scaled and histogramspecified CBCT images are compared to those in the corresponding CT image to determine how accurately the HU calibration of the images is recovered by matching the image histograms.
An automated analysis of the pixel values, described in Joshi et al (2017), is used to select and compare regions of pixels representing fat tissue in the CT and CBCT images. In brief: small circular regions of fat tissue in the CT image are automatically and randomly selected by identifying pixels with values between −100 and 0HU. The regions are then transferred to the CBCT image and are kept for analysis if they contain no air pixels and have a sufficiently small standard deviation, indicating that they contain pixels corresponding the one type of tissue. The mean value of all pixels in all regions of the image is then compared between CT and CBCT.
2.6. Image acquisition CT images used in this study were acquired using Siemens Somatom definition AS and Philips Brilliance Big Bore CT scanners. CBCT images were acquired using the Elekta XVI imaging system (XVI software version 4.5, Elekta, Crawley, UK). CT images of twelve patients, and three CBCT images per patient-36 CBCT images in total-were selected for analysis. All images were anonymised prior to processing and analysis.

Histogram peak identification
Example slices of the CT and CBCT pelvis images used to test the histogram transformation techniques are shown in figures 1(a) and (b) respectively 3 . Histograms created using pixels from all slices of the image volumes are shown in the lower panels.
The number and location of the peaks in the histogram, defined using equation (1), is susceptible to noise and depends on the bin width of the histogram. Figures 2(a) and (b) show the peaks identified in histograms of the CT image volume created using bin widths of 5HU and 30HU respectively.
Using a narrow bin width increases the smoothness of the histogram in areas where there are a significant number of pixels, for example the region around 1000HU where pixels represent soft tissue such as subcutaneous fat and muscle. However the noise in the histogram is also dramatically increased; there are far more peaks found which are not caused by any anatomical features. Figure 3 shows the final six steps of the procedure outlined in section 2.3 when applied to the histogram shown in figure 2(b), with n 2 target = . At each iteration of the procedure, the pair of peaks that are closest to each other, as defined by equation (2), are shown as solid red lines. The remaining peaks are shown as dashed blue lines. In the following iteration the two closest peaks are combined using equation (3).

Iterative peak combination
In the first thirteen iterations of the procedure the peaks corresponding to noise, located between 1750 and 4000HU, are combined to produce a peak at around 2350HU. It takes a further six iterations before the number of peaks have been reduced to n 2 target = , at which point we have a peak corresponding approximately to 'air' at around 50HU and one corresponding approximately to 'tissue' at around 1000HU. Figure 4(a) shows the results of applying the peak combination procedure to the histogram shown in figure 2(a). The procedure is robust and not affected by the large number of peaks produced by noise, other than needing a much larger number of iterations to combine the peaks until there are just  n 2 target = . The locations of the peaks obtained by algorithm are almost independent of the width of the bins used to create the histogram. For example the peak corresponding to tissue pixels has a HU value of 1009.5 when a histogram with bin width of 30 is used, and 1008.7 when a histogram of bin width 5 is used.
In figure 4(b) the iterative peak combination procedure is applied to the histogram of the CBCT image shown in figure 1(b). The air and tissue peaks are located after only a small number of iterations.

Histogram-transformed images
Having located air and tissue peaks on both the CT and CBCT image histograms (see figure 4), equation (4) can be used to define a linear transformation to map the CBCT image histogram peaks to those in the CT histogram. The same transformation is then applied to each pixel in the CBCT image volume. Figures 5(a) and (b) respectively show results of applying the iterative peak combination and histogram specification algorithms to the CBCT image shown in figure 1(b).  The histogram specification algorithm causes the resulting image histogram ( figure 5(b)) to very closely resemble the histogram used as a reference ( figure 1(a)). However, the histogram-specified image has reduced contrast and detail, and has additional artefacts compared to the histogram-scaled image shown in figure 5(a). Figure 6 shows the results of the pixel value analysis described in section 2.5. Each entry in the histogram represents the result of running the automated pixel value analysis on a particular image. The bin in which the entry is placed is determined by the average difference in fat tissue pixel values between the CBCT image and corresponding CT image. The pixels representing regions of fat in the raw CBCT images have pixel values that are on average around 136 HU different from the values in the corresponding CT images. When using the iterative peak combination algorithm to scale the histograms of the CBCT images the mean absolute difference in mean fat tissue pixel value is reduced to 23.7 HU. Note that after the iterative peak combination method has been used to identify relevant features of the histograms, the subsequent linear transformation applied to the CBCT image affects all pixel values equally. Therefore the pixels representing tissue other than fat, e.g. muscle, will be scaled in the same way.

Pixel value analysis
The histogram specification algorithm also produces CBCT images which have fat tissue pixels with values closer to those in the corresponding CT images.  With a mean absolute difference in mean fat tissue pixel value of 50.9HU the histogram-specified CBCT images do not resemble the CT images as closely as the CBCT images processed with the iterative peak combination algorithm. Additionally there are a small number of CBCT images which are poorer quality after being processed by the histogram specification algorithm, with mean differences in fat tissue pixel value of around 220HU. An example slice of one of these images is shown in figure 7. Despite having a histogram which closely matches the shape of the corresponding CT image histogram, the histogram-specified image ( figure 7(d)) exhibits many artefacts. Large areas of pixels representing fat tissue, along the posterior surface of the patient in particular, have values which are much too small after the application of histogram specification.
The results shown here for both the iterative peak combination and histogram specification methods were calculated using full CT and CBCT images without any initial registration or cropping. Differences in the imaged field of view (FOV) between CT and CBCT may contribute to the larger pixel value differences observed for the histogram-specified images, because the exact histogram shape will depend on the extent of the patient imaged. However the iterative peak combination method is not affected by differences in FOV because it only uses the positions of the major peaks in the histogram, not the full shape.

Discussion
In this report we present a novel method for analysing the histograms of CT and CBCT images, based on tools used to perform jet-finding in high energy particle physics. We demonstrated the use of this method to derive a scaling which matches CBCT to CT images based on their histograms. CBCT images that have been corrected in this way can be used to improve accuracy of radiotherapy dose calculations based on CBCT, or to improve image registration between CBCT and CT. We have used the iterative peak combination algorithm as an initial step in the shading correction of CBCT images (Marchant et al 2008). In this context the method has been shown to work robustly for images at different anatomical sites (pelvis, thorax, head and neck), and acquired with different fields of view .
CBCT images often contain artefacts and nonuniformities which result in pixels representing the same type of tissue having different values depending on their location in the image. The pixels, therefore, will occupy different bins of a histogram and will be affected differently by an algorithm such as histogram specification. The iterative peak combination procedure has been designed to keep in mind the relation between the features in CT and CBCT image histograms and the image structures which produced them.
When used in the context of jet-finding, equation (2) contains additional parameters, one of which controls the order that so-called pseudo-jets are combined: When p 1 = + we have our equation (2) and small peaks are combined before larger ones. Setting p 1 =would cause larger peaks to be combined before shorter ones, and with p = 0 we would ignore the height of the peaks and merge those that have HU values closest to each other. When used in jet-finding there are theoretically well-motivated reasons for the parameter p to take on each of the values, in different situations. We explored all three options in this work and have found that when iteratively combining the peaks of medical image histograms only a value of p 1 = + gives reasonable, robust results.
Our definition of a peak in a histogram (equation (1)) is sensitive to the bin width used to define the histogram, and thus the amount of noise present in the histogram, as illustrated in figure 2. However the iterative peak combination algorithm is robust against noise in the histograms, producing reliable results even when unrealistically large numbers of peaks are defined. The histogram specification algorithm is sensitive to noise in the histograms. A histogram of a real CT image is unlikely to have noise sufficient to obscure the peaks and important features of the histogram. However is it easy to imagine a case where the CT image histogram contains a large amount of noise that is transferred to the CBCT image histogram with undesirable results.
When defining our histograms we chose to use fixed-width bins. Variable-width bins could be used to provide more detail in regions of the histogram with higher bin frequencies, and reduce the impact of noise in regions with few entries per bin. The type of binning used would not affect the iterative peak combination procedure. Thing et al (2016) describe the use of histogram scaling as part of a procedure to recover HU accuracy in CBCT images. Smoothed histograms, created with 100 bins, were calculated using clinical and simulated CBCT projection images. Soft tissue and air peaks were identified and used to derive a calibration function similar to the linear function shown in equation (4). The number of bins used was sufficient to remove noise in the CBCT histograms and elucidate the important features-soft tissue and air peaks. However the same number of bins is unlikely to work with a CT image, and would potentially need to be tweaked if a CBCT projection image was simulated with a different exposure. The iterative peak combination procedure, in contrast, can be used to locate any number of interesting peaks and can be applied to histograms created with any number of bins. Arai et al (2017) use histogram specification to improve the accuracy of proton therapy dose calculations on CBCT images. The modified CBCT images resulted in improved proton dose calculation accuracy compared to raw CBCT. However the histogram-specified CBCT images of head and neck cancer patients appear to show reduced contrast between the fat and muscle tissue and other artefacts similar to those visible in figures 5(b) and 7(d). In addition, the histogram specification procedure was not used in regions where the CT and CBCT images were most different: at and below the level of the shoulders. Instead the CT and CBCT images were cropped to remove the problematic regions.
Park et al (2017) overcome a similar problem by applying histogram specification to each transverse slice of a CBCT image, after it has been registered to a corresponding CT image. The iterative peak combination procedure can also be used in situations such as this, by registering the CT and CBCT images and deriving linear transformations for each transverse slice. Figure 8 shows an example of a head and neck CBCT image before and after being processed with the iterative peak combination algorithm applied on a slice-by-slice basis. Severe artefacts in the region below the shoulders are reduced after processing. The slight exaggeration of the non-uniformity in some areas of the image is a consequence of the linear transformation applied to the image pixels, rather than the iterative peak combination method itself. A linear transformation derived by any means would produce the same effect, since the contrast of a CBCT image is typically less than that of a corresponding CT image. Therefore the CBCT histogram must be 'stretched' in order for its features to align with those in the CT image. This stretching also broadens the peaks in the CBCT image, which causes any non-uniformities to be exacerbated. However, after processing, the mean position of the soft tissue peak on each slice will be more closely matched to the CT image.
In addition to slice-by-slice processing, splitting the CBCT image into small three-dimensional patches might allow the histogram scaling approach to better correct the non-uniformities present in CBCT images, which cannot be removed using the global transformation defined in equation (4). Amit and Purdie (2015) use histogram specification to process CBCT images of patients undergoing radiotherapy treatment of the breast. They investigate both global and, in cases where relevant structures have been delineated on the CBCT images, ROI-based histogram specification. The corrected CBCT images are then used as part of a procedure which automatically generates radiotherapy treatment plans. The treatment plans generated were of a clinically acceptable quality. However the histogram-specified CBCT images produced using both the global and ROI-based correction exhibit quite severe artefacts. An approach such as the iterative peak combination algorithm is likely to produce images which perform equally well when used for radiotherapy dose calculations, but with fewer artefacts.

Conclusion
We present an iterative peak combination algorithm for analysing histograms of CT and CBCT images. The algorithm has only one parameter that determines the stopping point of the iterative combination and defines how many peaks are identified in the histogram. The algorithm has been used to automatically and robustly identify air and soft tissue peaks in CT and CBCT histograms. The peak locations are used to define a linear transformation that is applied to the CBCT image pixels. The resulting transformed image has pixels with values that are more closely calibrated to the HU scale. The proposed procedure has been compared to another commonly-used histogrambased transformation method: histogram specification. Iterative peak combination has been shown to be more robust that histogram specification and produces CBCT images with fewer artefacts and pixel values that more closely resemble those in corresponding CT images. Analysis of 36 CBCT images shows a mean difference in fat tissue pixel values between CT and raw CBCT of 136.7HU, compared to 50.9HU for histogram-specified CBCT images and 23.7 for those processed using the iterative peak combination algorithm. Figure 8. Iterative peak combination applied to each transverse slice of a head and neck image. (a) Shows a sagittal slice of a raw CBCT image of a head and neck cancer patient. Significant artefacts are present in the region of the image at and below the level of the patient's shoulders. The peak in the histogram at around 600 is caused by these pixels. (b) Shows the same image after each transverse slice is processed using the iterative peak combination algorithm. The HU values in the region around the shoulders have been brought in line with the remainder of the image, though some artefacts remain.