Calibration method for reducing ring artifacts in energy bin images for photon-counting ct

Abstract

An apparatus and method for calibrating a photon-counting computed tomography (PCCT) apparatus, includes performing for each offset from a rotation center, a scan of a calibration phantom placed at the offset to obtain calibration data including measured count rate data. Location parameters of the calibration phantom are estimated for each of the plurality of offsets, based on the obtained calibration data. Thickness data of the calibration phantom is determined based on the estimated location parameters. A curve is estimated, by curve fitting using the measured count rate data and the determined thickness data, indicating a relationship between the energy bin count rate and phantom thickness, for each energy bin of each pixel of the plurality of pixels. A first look-up table is generated for pixel-to-pixel count rate variation correction, and the first look-up table is applied to acquired energy bin count data to reduce ring artifacts in an image reconstructed using the acquired energy bin count data. A second look-up table is generated for both pixel-to-pixel count rate variation correction and additional count-rate non-linearity correction, using water as the calibration material, and the second look-up table is applied to reduce both ring artifacts and cupping or doming artifacts in an image reconstructed using the acquired energy bin count data.

Description

[0001] CALIBRATION METHOD FOR REDUCING RING ARTIFACTS IN ENERGY BIN IMAGES FOR PHOTON-COUNTING CT

[0002] CROSS-REFERENCE TO RELATED APPLICATIONS

[0003] This application claims the benefit of priority to provisional Application No.

[0004] 63 / 726,452, filed November 29, 2024, the entire contents of which are incorporated herein by reference.

[0005] BACKGROUND TECHNICAL FIELD

[0006] The present disclosure is directed to a calibration framework for photon counting detectors in computed tomography imaging. In particular, the calibration framework involves the generation of two look-up tables: one for pixel-to-pixel variation correction and the other for both pixel-to-pixel variation and additional count-rate non-linearity correction based on water.

[0007] DESCRIPTION OF RELATED ART

[0008] Medical imaging involves obtaining images of the body. A computerized tomography (CT) scan is a type of imaging that uses X-ray techniques to create detailed images of the body, creating cross-sectional images, also called slices, of the bones, blood vessels, and soft tissues inside the body.

[0009] In computed tomography (CT) imaging, photon-counting detectors (PCDs) offer many advantages over conventional energy-integrating detectors (EIDs). PCDs directly convert X-ray photons into electrical signals, counting individual photons and their energy levels, which leads to better image quality and diagnostic capabilities compared to EIDs. A

[0010] PCD classifies incident photons into several energy bins based on their energy by comparing all pulses against several different thresholds. However, PCDs present challenges, one of which is pixel-to-pixel variation in detector outputs. Individual pixels in PCDs exhibit greater variability compared to those in EIDs, leading to more pronounced ring artifacts in CT images. Moreover, variability between pixels can depend on the energy bin of the PCD. Also, spectral distortion impacts each bin differently, making correction more difficult. As such calibration is a complex process for PCDs.

[0011] Ring artifacts in photon-counting detectors are concentric circles in reconstructed images caused by mis-calibrated or defective pixels on the detector array. These artifacts appear as rings because a misbehaving pixel will produce a consistently erroneous reading at each angular position as the detector rotates around the object. This can lead to a degradation of image quality, potentially affecting diagnostic accuracy.

[0012] Moreover, the variability between pixels differs across energy bins of PCDs, which further complicates the issue. In PCDs, the variability between pixels differs across energy bins primarily due to two physical effects, referred to as charge sharing and pulse pile-up. These effects, which are influenced by the size of the pixel and the energy of the incident photons, cause the measured count rate and energy to vary from pixel to pixel in an energydependent manner. To maximize the benefits of PCDs, it is crucial to address pixel-to-pixel variations and minimize ring artifacts.

[0013] Efforts to address pixel-to-pixel variations have pursued both hardware and software approaches. Hardware efforts include offset correction techniques, improved crystal technology, and advanced application-specific integrated circuits (ASICs) design. Hardwarebased solutions provide more intrinsic, long-term improvements, but they tend to be more complex, costly, and time-consuming to implement. Although hardware advancements are essential and should continue, fully eliminating pixel-to-pixel variations is unlikely in the near future. Therefore, software-based methods should remain and complement hardware solutions. Prior to the advent of PCDs, various correction methods using EIDs were extensively studied. These methods can be broadly categorized into three groups: sinogram -based methods, image-based methods, and calibration methods. Both sinogram-based and imagebased approaches generally follow similar steps, first identifying lines in the sinogram or image domain (after polar coordinate transformation) using different algorithms, such as median or mean filters, Fourier transform -based algorithms, wavelet transform-based algorithms, gradient-based algorithms, and more recently, deep learning-based algorithms. After detecting the artifacts, these methods apply correction. While these methods are simple and flexible, their effectiveness can significantly vary depending on the algorithm's ability to accurately detect and correct artifacts based on local patterns. Some of these methods that lack precision in detecting and correcting artifacts often compromise image resolution.

[0014] Various calibration methods address pixel-to-pixel variations by empirically estimating the relationship between pixel outputs and influencing parameters, such as tube current and calibration material thickness, on a per-pixel basis. These methods have often been referred to as “flat-field correction” in previous studies, although they address both pixel-to-pixel variations and beam intensity normalization. More integrated calibration methods also aim for other objectives, such as beam-hardening correction (BHC) or material decomposition (MD), which involve calibration with one or more materials and incorporate pixel-to-pixel variation correction as part of the overall process. These calibration methods can achieve highly accurate corrections without compromising image resolution when calibration data includes comprehensive coverage of X-ray intensities and spectra. Acquiring such extensive calibration data for all potential coverage of X-ray intensities and spectra, however, is very challenging. Therefore, it is important to balance calibration efficiency with calibration performance for pixel-to-pixel variation correction. For EIDs, a moderate level of calibration data generally suffices for pixel -to-pixel variation correction, as calibration using two or multiple condition points for each pixel is often sufficient, and the calibration process remains relatively simple. In contrast, PCDs, where they exhibit higher pixel-to-pixel variation and non-linear responses influenced by factors like charge sharing and pulse pileup, require denser calibration data to account for added complexity, alongside a more sophisticated calibration process.

[0015] An object of the present disclosure is a calibration framework tailored for PCDs, that addresses the challenges posed by their increasingly complex behavior. A further object is to maximize the benefits of PCDs, by addressing pixel-to-pixel variations and minimizing ring artifacts. A further object is a streamlined calibration data acquisition strategy that uses CT scans of calibration phantoms combined with a robust calibration process.

[0016] Aspects of this technology are described in the article: D. Lee, X. Zhan, W. Yang Tai, S. Subramanian, W. Zbijewski and K. Taguchi, " Calibration Method for Reducing Ring Artifacts in Energy Bin Images for Photon-Counting CT," in IEEE Transactions on Medical Imaging (2025), which is incorporated herein by reference in its entirety.

[0017] Aspects of this technology are described in the article: Donghyeon Lee, Xiaohui Zhan, W. Yang Tai, Wojciech Zbijewski, Katsuyuki Taguchi, " Pixel-to-pixel variation correction using cylindrical phantoms in photon-counting CT: total count results," Proc. SPIE 13405, Medical Imaging 2025, which is incorporated herein by reference in its entirety.

[0018] SUMMARY

[0019] An exemplary embodiment is a calibration method that includes both a comprehensive data acquisition method using CT scans and a three-step calibration process to address ring artifacts in reconstructed energy bin images for photon-counting CT. The data acquisition method is a streamlined data acquisition strategy using CT scans of calibration phantoms, enabling comprehensive phantom thickness data acquisition across all pixels. The three dedicated calibration steps generate two look-up tables: one for pixel-to-pixel variation correction and the other for both pixel-to-pixel variation and additional count-rate nonlinearity correction based on water. The calibration method serves as an effective solution for addressing pixel-to-pixel variation in PCDs.

[0020] An aspect is a method for calibrating a photon-counting computed tomography (PCCT) apparatus. The method can include performing, using the PCCT apparatus, for each offset position of a plurality of offset positions from a rotation center, a scan of a calibration phantom placed at the offset position to obtain calibration data including measured count rate data, wherein the PCCT apparatus is configured with a semiconductor-based photon-counting detector (PCD) having a plurality of pixels; estimating, by processing circuitry, for each of the plurality of offset positions, location parameters of the calibration phantom based on the obtained calibration data; determining, by the processing circuitry based on the estimated location parameters, thickness data of the calibration phantom; estimating, by curve fitting using the measured count rate data and the determined thickness data, a curve indicating a relationship between the energy bin count rate and phantom thickness, for each energy bin of each pixel of the plurality of pixels; generating, by the processing circuitry based on the estimated curves of each energy bin of each pixel, a first look-up table for pixel-to-pixel count rate variation correction; and applying, by the processing circuitry, the first look-up table to acquired energy bin count data to minimize artifacts in an image reconstructed using the acquired energy bin count data.

[0021] A further aspect is an apparatus for calibrating a photon-counting computed tomography (PCCT) device. The apparatus can include processing circuitry configured to cause the PCCT apparatus to perform, for each offset position of a plurality of offset positions from a rotation center, a scan of a calibration phantom placed at the offset position to obtain calibration data including measured count rate data, wherein the PCCT apparatus is configured with a semiconductor-based photon-counting detector (PCD) having a plurality of pixels, estimate, for each of the plurality of offset positions, location parameters of the calibration phantom based on the obtained calibration data; determine based on the estimated location parameters, thickness data of the calibration phantom; estimate, by curve fitting using the measured count rate data and the determined thickness data, a curve indicating a relationship between the energy bin count rate and phantom thickness, for each energy bin of each pixel of the plurality of pixels; generate, based on the estimated curves of each energy bin of each pixel, a first look-up table for pixel-to-pixel count rate variation correction; apply the first look-up table to acquired energy bin count data to minimize artifacts in an image reconstructed using the acquired energy bin count data.

[0022] The foregoing general description of the illustrative embodiments and the following detailed description thereof are merely exemplary aspects of the teachings of this disclosure, and are not restrictive.

[0023] BRIEF DESCRIPTION OF THE DRAWINGS

[0024] A more complete appreciation of this disclosure and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:

[0025] FIG. 1 is a flowchart of a calibration method, in accordance with exemplary embodiments;

[0026] FIG. 2A is an illustration of a CT system scanning a phantom with offsets;

[0027] FIG. 2B illustrates a graph of water thickness sampling distributions from a single CT scan;

[0028] FIGs. 3A, 3B, and 3C illustrate different shapes of phantoms; FIG. 4A illustrates template images;

[0029] FIG. 4B illustrates a weight image;

[0030] FIG. 4C is a scatter plot of water thickness versus count-rate obtained by data acquisition on a fitted curve;

[0031] FIG. 4D is detector row data with float water thickness of 10 cm;

[0032] FIG. 5A is a graph of water thickness to count-rate for three neighboring pixels;

[0033] FIG. 5B is a graph of correction coefficients with respect to the line integral for three pixels;

[0034] FIG. 6A illustrates a head phantom;

[0035] FIG. 6B illustrates a multi-energy body phantom;

[0036] FIG. 7 illustrates an example of a photon-counting Computed Tomography (PCCT) scanner system;

[0037] FIG. 8A is a plot of optimal pose parameters;

[0038] FIG. 8B is a graph of collected data from the optimal pose parameter;

[0039] FIG. 9 illustrates graphs of curve fitting results of energy bin count-rates against water thickness for a peripheral pixel and center pixel;

[0040] FIG. 10 illustrates graphs of good pixel fitting results;

[0041] FIG. 11 illustrates reconstructed images of a 24 cm water phantom without offset, obtained using different calibration methods;

[0042] FIG. 12 illustrates reconstructed images of a 24 cm water phantom with a 5 cm offset, obtained using different calibration methods;

[0043] FIG. 13 illustrate graphs of radial mean profiles of total count, bin 1, and bin 5;

[0044] FIG. 14 illustrate graphs of high-frequency and low-frequency RMSD of total count and all energy bin counts; FIG. 15 illustrates reconstructed images of a head phantom, obtained using different calibration methods;

[0045] FIG. 16 illustrates reconstructed images of a multi-energy body phantom, obtained using different calibration methods;

[0046] FIG. 17A is a graph of average ROI differences; and

[0047] FIG. 17B illustrates zero-count occurrences in the total count and bin 2.

[0048] DETAILED DESCRIPTION

[0049] In the drawings, like reference numerals designate identical or corresponding parts throughout the several views. Further, as used herein, the words “a,” “an” and the like generally carry a meaning of “one or more,” unless stated otherwise.

[0050] Furthermore, the terms “approximately,” “approximate,” “about,” and similar terms generally refer to ranges that include the identified value within a margin of 20%, 10%, or preferably 5%, and any values therebetween.

[0051] An aspect of the present disclosure is a photon-counting detector (PCDs) in computed tomography (CT) imaging that addresses a pixel-to-pixel variation in order to minimize ring artifacts in CT images. Moreover, disclosed is a framework that corrects variability between pixels depending on the particular energy bin of the PCD, and thereby reduces the impact of spectral distortion for each bin. The calibration framework leverages gantry-rotating CT scans of calibration phantoms placed at multiple offsets from the rotation center, enabling a denser sampling of phantom thickness data across all pixels. The framework includes a streamlined calibration data acquisition strategy that uses CT scans of calibration phantoms, combined with a calibration process to achieve higher calibration accuracy.

[0052] FIG. 1 is a flowchart for a calibration method in accordance with exemplary embodiments. The disclosed calibration method includes two primary stages: step S102, acquisition of calibration data, and step S104, a calibration process. The three steps of the disclosed calibration process, S104, include a step S104a of phantom location estimation, a step S104b of curve fitting, and a step S104c of look-up table generation. The results of the calibration process S104 can be applied in an application step S106.

[0053] Calibration data acquisition step (S102)

[0054] The data acquisition method includes positioning calibration phantoms with multiple offsets from the iso-center, and conducting CT scans. The calibration data acquisition using the plural offsets is a technique that allows for the acquisition of comprehensive phantom thickness data for the entire detector pixels using CT scans. One method of calibration data acquisition is described in Lee et al., “X-Ray Transmittance Modeling-Based Material Decomposition Using a Photon-Counting Detector CT System,” IEEE Trans. Radiat. Plasma Med. Sci., vol. 5, no. 4, pp. 508-516, 2021, incorporated herein by reference.

[0055] FIG. 2A and FIG. 2B illustrate an off-centered cylindrical phantom in a CT system and water thickness sampling distributions from a single CT scan of a 24-cm-diameter cylindrical water phantom at different offsets, respectively. Depending on the offsets, the sampling range with respect to water thickness varies significantly. As the offset increases, the sampling area expands in the channel direction while the sampling density in the water thickness direction decreases. The sampling area continues to expand until the center of the cylindrical phantom reaches the boundary of the field-of-view (FOV). Note that due to the cylindrical shape, the sampling density is highest around the diameter of the water phantom and decreases as the thickness reduces.

[0056] In one embodiment, data is collected using three different-sized water phantoms at three different offsets, resulting in a total of nine datasets, to achieve sufficient coverage and sampling density with respect to water thickness across the entire detector array. In this embodiment, water is used as the calibration material.

[0057] In one embodiment, the shape of a phantom is a cylinder, as in FIG. 3 A. Other phantom shapes can include an ellipsoidal cylinder, as in FIG. 3B, and a cone, as in FIG. 3C.

[0058] In one preferred embodiment, the calibration phantom is water because water phantoms have similar attenuation properties to human tissue. The phantom materials can vary depending on target objects. Other phantom materials can include plastics, such as acrylic, polystyrene, polyethylene, aluminum, and bone.

[0059] Calibration process (S104)

[0060] The calibration process includes the following three steps.

[0061] Step S104a is a phantom location estimation step that estimates the phantom location parameters. These parameters are used to accurately calculate the phantom thicknesses corresponding to measured count-rates.

[0062] Step S104b is a curve fitting step that uses curve fitting to estimate the relationship between measured count-rate and phantom thickness. The resulting curve defines the expected outputs of each pixel with respect to phantom thickness. Achieving accurate curve fitting is needed to ensure reliable corrections.

[0063] Step S104c is a look-up table generation step that generates look-up tables (LUTs) for correction coefficients. Two types of LUTs are generated: the first corrects only for pixel -to-pixel variation, while the second addresses both pixel-to-pixel variation and count-rate nonlinearity with respect to water thickness.

[0064] The calibration process involves two types of data: total counts and energy bin counts. While all steps are conceptually the same with either type of data, total count calibration uses total count-rates, whereas energy bin count calibration uses count-rates of each energy bin for Steps S104b and S104c. Step S104a is identical for both types of data as they both use total counts for estimation. For the sake of brevity, the following steps are described in the context of energy-bin count calibration.

[0065] Phantom Location Estimation, SI 04a

[0066] This step determines optimal location parameters and calibration phantom thicknesses using alignment of 3D template projections with line integrals. A template image is generated by creating a 3D image of an object to be scanned (i.e., calibration phantom). The 3D image is created using information about the CT scanner’s geometry and the shape of the object when it is perfectly aligned (without any shifting). The template acts as a reference to help match the template with the real measurement that is collected from the scanner.

[0067] The location parameters for the x and y axes in the image domain are represented as a vector r = (rx, ryThe optimal parameters align projections of the calibration phantom’s 3D

[0068] template image pswith line integrals ls. The 3D template image is generated using prior knowledge of the CT system geometry and the calibration phantom. FIG. 4A shows sample example slices of the 3D template images. The line integral is calculated as

[0069] , »..m ax f y *f / a & Ix..

[0070] S = -in ( - ’X t- * + ( 1 )

[0071]

[0072] where y^Ui Vi 0i sand y

[0073]

[0074] “"ti vare the measured total count-rate and air count-rate, respectively,

[0075] y

[0076]

[0077] *, u, v,e,s = Sb yb, u, v, e, s and y“^v= Sb yfiu,v■ yb, u, v, e, s and yf\vare the measured energy bin count-rate and bin air count-rate, where 6 G {1, 2,...,bmaxu G (1, 2, umax]

[0078] , v £ {1, 2,..., vmax), 0 E {1, 2,..., 0max) and s G {1, 2,.... smax} denote the energy bin, pixel row, pixel column, projection angle, and dataset indices, respectively. The max function chooses the larger value between two inputs, e and 8 are small positive values added to prevent a logarithmic singularity when the ratio is close to zero and to reduce zero count bias from logarithmic transformation, respectively. A weight image is generated to help improve the alignment process. This sub-step focuses on certain parts of the template image that are most important. A weight image is created that gives higher importance to important regions, such as edges, and less importance to other parts. This sub-step helps to minimize the impact of unimportant areas, ensuring that the final matching is precise.

[0079] The weight image is employed to isolate edge volumes that are highly related to the location parameters and to reduce the nonlinearity between the projections of the 3D template image and the line integrals. The weights are defined as follows:

[0080] _ f 1 iffmin —rmax / .-JA

[0081]

[0082] - |;() Qtherwlse’

[0083] where Tminand Tmaxare the minimum and maximum thresholds, respectively. An example of the weight image generation is presented in FIG. 4B.

[0084] To estimate the optimal location parameters, an iterative algorithm uses a weighted normalized cross correlation (NCC) cost function. The NCC cost function is described in J. P. Lewis, “Fast Normalized Cross-Correlation,” Vis. Interface, 1995, incorporated herein by reference in its entirety. The use of the NCC is motivated by two main reasons: normality, ensuring that the similarity measure is not affected by absolute intensity differences, and pixel-by-pixel comparison, leveraging spatial information for precise alignment.

[0085] The goal of estimating the optimal location includes finding the best alignment between the reference template and the real scan data by adjusting the position of the template. An optimization method iteratively estimates the offset parameters of the template image until it matches the measured data as closely as possible. The quality of the match is measured using the NCC cost function, which determines how similar the two images are, ignoring any overall brightness differences. Once the best fit is found, the corresponding alignment information is used to estimate the thickness of the object for further calibration steps.

[0086] The weighted NCC, which serves as the similarity measure for the cost function, is defined as:

[0087] lVCCw(a, b. ) = / sV f J fcW; » fed s,-fJ.— ajt&jf..; A I.,{3)

[0088]

[0089] \ ’ J where a and b are the weighted means of the 3D matrices a and b: a = — ~ — — and b =

[0090]

[0091]

[0092] W

[0093] i] kW k• The optimal location parameters tor each dataset s are obtained by solving:

[0094] rs* — argmin CC w.8;dS£(Tr(ps)Y (4) where L is the forward projector generating simulated line integrals based on system geometry, and Tris the 2D translation function defined by the location vector r.

[0095] Provided the optimal location vector r*, the calibration phantom thicknesses (mm) can be calculated for each dataset:

[0096] t

[0097]

[0098] s= £(?,,. (ps)), (5) Then, the determined phantom organizes thicknesses and the measured count-rates as follows:

[0099] (t£h. U - U t„, (6)

[0100]

[0101] U-Uwhere the sort function merges the third and fourth dimensions of a matrix and arranges the merged elements in ascending order. The sort (i) also merges input matrix dimensions but

[0102] arranges data based on the sorted order of t^.. n G {1,2,..., nmax] is a new data index.

[0103] Curve Fitting, SI 04b It has been found that Gaussian averaging alone does not result in a sufficient curve. A curve fitting step is applied to obtain improved straightness. The curve fitting step consists of three sub-steps for curve fitting: (a) Gaussian-weighted averaging, (b) curve fitting, and (c) interpolation. The curve fitting step is conducted independently for each energy bin within each pixel. The indices b, u, and v are omitted for simplicity:

[0104]

[0105] and The measured data points are smoothed using a method that gives more weight to nearby points (Gaussian-weighted averaging). This helps reduce noise in the data without distorting the values too much. A Gaussian-weighted averaging is applied with a small sigma to standardize the non-uniformly distributed data with respect to water thickness:

[0106] Y, frf1). (1) \ (i)—Z; i) fi) 7 ■> (1)

[0107] w

[0108]

[0109] here wgauss(a, b; a) = ^=e-Ca~b)2 / 2rT2and e (t^, t^,... J.m et1’2- -mmax] is a new data index where data density with respect to water thickness is uniform like A

[0110]

[0111] fC) = t^+1— for all m.

[0112] Curve fitting can be performed using locally weighted scatterplot smoothing (LOWES S). The LOWES S process fits a local second-order polynomial using weighted data, emphasizing nearby points, to produce a smooth curve. A basic version of LOWES S is described in W. S. Cleveland, “Robust Locally Weighted Regression and Smoothing Scatterplots,” J. Am. Stat. Assoc., vol. 74, no. 368, pp. 829-836, 1979, incorporated herein by reference in its entirety. This technique adjusts how much data it uses depending on the shape of the curve. Given two vectors, a and b, along with a span a, the LOWESS function is defined as:

[0113]

[0114] = LOWESS(aifa, b; a), (8) where a and b are one dimensional vectors: a = (a, and b = (b,...bimnx). The span a is a parameter ranging from 0 to 1 that determines the proportion of neighboring data points to the entire data points used for smoothing. LOWESS is further described below.

[0115] To enhance LOWESS fitting, two strategies include logarithmic transformation and dynamic span. The logarithmic transformation simplifies the fitting task by making the relationship between count-rate and thickness closer to linear. The dynamic span is used to account for the variation in count-rate properties based on phantom thickness (refer to FIG.

[0116] 4C). For a thin thickness, the curvature of the logarithmic count-rates is higher, indicating sharper changes, with lower noise, making a smaller span more appropriate. For a thick thickness range, the curvature is lower, indicating smoother changes, with higher noise levels, making a larger span more suitable. The curve fitting can be expressed as follows:

[0117] = exp t(1), a^yn)], (9)

[0118] where ft1) G ( x.t^ i,...,, Lmax )'' and Iny17W E ( vIny ✓ ^ x, Iny s ^ x.,... Iny y ^ "hnax' ). a^,, c nis a

[0119] dynamic span defined as:

[0120] d m T" ^min» (10)

[0121]

[0122] where aminand amaxare the minimum and maximum span values, respectively, tp denotes

[0123] the boundary thickness (mm) beyond which the maximum span is applied.

[0124] Once there is a smooth curve, it is used to estimate the count rates at evenly spaced thickness points. In this way, the results are easier to compare and analyze.

[0125] In particular, based on the outputs of Eq. (9), the fitted count-rates are calculated at uniformly spaced thickness intervals using interpolation as:

[0126] fg'j f l';

[0127] (2) — (1).—z — (1) — (l)x ♦ < •* Jjg—(i) ( -0

[0128]

[0129]

[0130] where = (q — 1) X At for < t^+1■ q E {1, 2,..., qmax} is a new data index and At is the thickness interval (mm). For t_qA((2)) that exceeds the maximum thickness (i.e., t(2)> t(1))

[0131] If the thickness exceeds the measured range, the curve is extended to guess the count rates for those thicknesses (extrapolation). Extrapolation is applied using a fourth order polynomial.

[0132] Look-up Table Generation, SI 04c

[0133] The calibration method next includes generating each of the two LUTs individually. Even when different pixels are measuring the same spot in the object, they show slightly different count-rates because of response variations. To solve this, the method creates a look-up table that corrects each pixel’s measurement to match the behavior of an “good” pixel. A good pixel here refers to an “ideal” pixel that exhibits no pixel -to-pixel variation, with outputs depending solely on the incident X-ray intensity and spectrum, while still having physical effects, such as beam hardening, scattering, charge sharing, and pulse pileup.

[0134] Pixel-to-pixel correction maps measure count-rates to ’’good pixel” count-rates. The behavior of the good pixel is estimated by fitting a smooth curve of the pixel outputs at a specific phantom thickness. In Step 104b, the curve estimated for each pixel represents its expected response to water thickness. A set of these curves shows a variation, as shown by the curve 402 in FIG. 4D. If all pixels were good pixels, the count-rate profile of the pixels would form a smooth bell curve like the curve 404 in FIG. 4D, as the incident x-ray intensities and spectra reaching the pixels vary smoothly with pixel position. Estimating the good pixel behavior basically involves identifying this smooth bell curve across water thicknesses. The X in FIG. 4D at 406 marks a point that corresponds to a point in FIG. 4C, marked by X at 406. (2')

[0135] First, moving averages of y^u,v qare calculated based on a 2D window. The size of

[0136] the 2D window is defined based on the detector block size, as each pixel block may behave similarly. The calculation is as follows:

[0137] „IH-WF / 3 -, M+W'a / 2 (3)

[0138]

[0139] where Wuand Wvrepresent the horizontal and vertical sizes of the 2D window, respectively. Then, the smooth curve is calculated using the LOWESS function as the following equation:

[0140]

[0141] = “P u. in '■ “»))’ where l

[0142]

[0143] nyg°v q= (Inyg^v q,....\ny^Umax,v, <?)• U is a horizontal pixel index vector: U

[0144] = (1, umax). agis a span for the good pixel estimation.

[0145] The coefficient for pixel -to-pixel variations is calculated as follows:

[0146] PPF / (2)n4x

[0147]

[0148] The correction for water linearization makes sure that when measuring different thicknesses of water (which has a predictable response to X-rays), all the pixels show the expected results. Water is used as a reference because its response to X-rays is well understood and consistent. This correction not only addresses the differences between pixels, but also ensures that the response is linear with respect to water thickness. It makes all pixels behave in the same way when measuring water, which is crucial for accurate imaging.

[0149] Water linearization maps measure count-rates to count-rates whose line integrals are linear with respect to water thickness. Note that this correction includes both water linearization and pixel-to-pixel variation correction as it ensures that all pixels behave in the same way. The energy bin count-rate with water linearization is calculated as follows:

[0150] J

[0151]

[0152] JX =xexpC-^fEfc) X (<? - l)dt), (15) where iwEb) is a water attenuation coefficient (cm-1) at the x-ray energy Eb■ keV). Ebis the selected x-ray energy to represent energy bin b.

[0153] The coefficient for water linearization is calculated as follows:

[0154] W=(3&), (2)

[0155]

[0156] J b,u,v,q f b,u,v,q ’ (16)

[0157] The result is two look-up tables: one for fixing pixel-to-pixel variations and another for making sure water thickness measurements are consistent. Using these tables helps produce high-quality images where all parts of the image show a uniform response, free from unwanted pixel-to-pixel differences.

[0158] Application (SI 06)

[0159] A test method is used to apply the look-up tables (LUTs) to the test data. A goal of the test method is to find where the test data falls between two reference count rates and use linear interpolation to determine the correction coefficient.

[0160] C2')

[0161] The count-rate ybiu,v,q obtained from the fitting in Step 2 of the calibration steps can

[0162] be applied to the correction coefficients in the LUTs to test data

[0163]

[0164] v 0.

[0165] First, the process identifies the index in which a line integral of the test data count-rate

[0166] Z^v efalls between two consecutive reference line integral Z),vand Z),r,q+i, such that

[0167] v.q G V.e < G v.q+1- Here,thereference line integrals are defined as zg^v q= -

[0168] In (yS,v,Q / y£,v,i + 0 andlb<, l,v,q+i= - In (yS,v,Q+i / y£,v,i + 0- Then, linear interpolation of the coefficients from Eq. (14) is performed to compute the pixel-to-pixel

[0169] variation-corrected count-rate y^u.’v,^6 orEq. (16) to compute the water-linearized count¬

[0170] rate yb!?v ’ respectively. Once the look-up tables are generated, they can be evaluated using test data. The testing method can use linear interpolation to determine the correction coefficient.

[0171] First, an interpolation ratio is calculated given a measured line integral for test data as

[0172] pest _ [LUT

[0173] _ 1

[0174] / L —jLUT > iLUT

[0175]

[0176] lq+llq

[0177] where ltestis the measured test line integral. lqUT< ltest<

[0178]

[0179] q e {1, 2,..., qmax{

[0180] is an index for look-up tables.

[0181] Then, the pixel-to-pixel variation-corrected (PPVC) count rate and water linearized (WL) count rate are computed using:

[0182] ,test, PPVC > frPPVC I i (rPPVC „PPVc\l..test

[0183] —cq -t- * — CqJj * y

[0184] and

[0185] yjtest. WL — i 3 # (rWL..test y

[0186]

[0187] - icq+ 7v?+i~cq j) y

[0188] where ytestis the measured test count rate. cppycand c^Lare the coefficients from

[0189] the LUTs for pixel-to-pixel variation correction and water linearization, respectfully.

[0190] The test method may alternatively be perform according to one of three variations. In a first variation, the test method is performed as cascaded LUTs. In this variation, the interpretation ratio for the application of the calibration method can be calculated by either the measured line integral or the measured count rate as follows:

[0191] (1)=

[0192] tOf

[0193] q4-l

[0194] (2) A =

[0195]

[0196] The performances of the two approaches are comparable and similar. Depending on how the look-up tables are created, both are available for application use cases.

[0197] In a second variation, the test method is conducted sequentially. In this variation, the pixel -to-pixel correction and water linearization can be conducted sequentially as:

[0198] ."test. PPVC — frPPV I i * (rPPVC PPVcAl..test y

[0199]

[0200] — \cq+ A* \C<7+1 ~Cq J) V and

[0201] ytest, WL=^CWL+^PPVC * >CWL) } * ytest. PPVC

[0202]

[0203] where, LUT — • For this sequential approach, the look-up table for water

[0204]

[0205] linearization is built to be applicable for the sequential correction. As in the first variations, for the sequential approach, the input for the interpolation ratio calculation can also be changed for both the measured line integral of the measured count rate.

[0206] In a third variation, direct conversion to line integrals can be performed as:

[0207] ltestlPPVC=^cPqPVC+^CPPVC _CPPVC)} *ytest

[0208]

[0209] and

[0210] Itest. WL= +X * *ytest

[0211]

[0212] The third variation has an advantage in computational speed compared to the other variations.

[0213] FIGs. 5A and 5B illustrate the differences between two correction results and coefficient values. FIG. 5A shows water thickness-count-rate plots for three neighboring example pixels 502, 504, and 506, a good pixel 508, and a water-linearized (WL) 510 pixel. The good pixel curve (508) passes through the average positions of the neighboring pixels, while the WL pixel line (510) maintains linearity with respect to water thickness, showing their distinct purposes.

[0214] FIG. 5B presents the corresponding correction coefficients as a function of line integral.

[0215] Applications of calibration methods with 2D interpolation that use different mAs. Advanced CT systems have a feature that modulates tube current while performing CT scans, called Tube Current Modulation (TCM), to reduce exposure dose to patients. However, this introduces a challenge of acquiring calibration data and constructing look-up tables (LUTs) for every possible mA value, which is impractical due to the continuous variation of mA during scanning.

[0216] For mA values that are not explicitly included in the LUTs, this limitation is addressed by applying two-dimensional interpolation. This approach simultaneously considers both count rates (or line integrals) and tube currents (mA) to interpolate correction coefficients. Specifically, the LUTs are constructed for discrete mA values and specific count-rate intervals, and for any given test values, the nearest neighboring mAs and countrates in the LUT are identified. Next, an interpolation method is used to compute the correction coefficients by combining the four neighboring LUT coefficients corresponding to these test values. The interpolation method can include, but is not limited to bilinear, bicubic, or B-spline interpolation.

[0217] Multi-point Flat Field Calibration Method

[0218] A multi-point flat field calibration method is provided as a comparison method. This calibration approach involves acquiring open field data, unobstructed X-ray images, at different X-ray intensities. The X-ray intensity is controlled by varying the X-ray source filament current (mAs levels). The calibration method using open field data using EIDs and PCDs is described in Seibert et al., Kwan et al., Schmidgunst et al., Altunbas et al., and Lifton et al. These works demonstrate its effectiveness in addressing pixel-to-pixel variations, making it a suitable comparison to the disclosed method.

[0219] The calibration process consists of the following steps:

[0220] Step 1 (Acquisition of open field data)'. Acquire open field images at different X-ray tube currents by adjusting the source's mAs. Then, calculate the mean intensity across all pixels at each x-ray tube current.

[0221] Step 2 (Curve fitting) '. For each pixel, plot the pixel intensities against the mean intensities. Use piecewise linear or polynomial fitting to model the relationship, providing a correction coefficient for each intensity range.

[0222] Step 3 (Generating look-up tables) '. Generate LUTs based on the fitting results, storing correction coefficients.

[0223] In Step 3, a piecewise linear interpolation is applied to calculate the correction coefficient within the calibration data points and second-order polynomial extrapolation for points outside the calibration range.

[0224] Experimental Assessment

[0225] Next, an outline of the experimental setup and the quantitative evaluation metrics are described.

[0226] (1) Experimental setup

[0227] Imaged objects'. An experimental setup uses three different-sized water cylindrical phantoms (18, 24, and 32 cm in diameter), a head phantom, and a multi-energy body phantom (40 cm x 30 cm). The head phantom is made of the same resin-based material, but different densities (see FIG. 6A). The ellipsoidal body phantom consists of a base material (solid water) and 12 inserts with varying iodine and calcium concentrations, as well as blood, brain, and adipose tissue mimicking materials (see FIG. 6B).

[0228] Data acquisition conditions'. An X-ray tube was operated at a tube voltage of 120 kVp with tube currents 10, 50, 100, and 200 mA, and a bowtie filter was installed in front of the X-ray tube. The three differently sized water phantom data were scanned with 10, 15, and 18 cm offsets at tube current 200 mA and used for calibration, while 24 cm water phantom data were scanned with 0 cm and 5 cm offsets at tube current 200 mA and used for testing. Both the head and the body phantoms were tested. In each CT scan, we acquired 1200 projections performed over one second.

[0229] PCD'. The experimental setup includes a prototype PCD-based CT (Canon Medical Systems, USA). See X. Zhan, R. Zhang, X. Niu, I. Hein, B. Budden, S. Wu, N. Markov, C. Clarke, Y. Qiang, H. Taguchi, K. Nomura, Y. Muramatsu, Z. Yu, T. Kobayashi, R.

[0230] Thompson, H. Miyazaki, and H. Nakai, “Comprehensive evaluations of a prototype full field-of-view photon counting CT system through phantom studies,” Phys. Med. Biol., vol. 68, no.

[0231] 17, 2023, incorporated herein by reference in its entirety. The detector is a curved detector, which consists of 2688x192 micro pixels, and uses cadmium zinc telluride as the crystal material. Threshold energies are set at 30, 45, 55, 65, and 80 keV. Anti-scatter grids are placed on the detector.

[0232] Calibration parameter settings'. The parameter settings used in the experimental setup are as follows: in Step 1, a = 0.1 mm and At^^O.l mm for Gaussian weighted moving averaging; in Step 2, amin= 0.1 mm for Gaussian weighted moving averaging; in Step 2,amin=005, amax= 0.8, and tB= 320 mm for the dynamic span configuration, and At = 0.1 mm for interpolation; and in Step 3, ag= 0.3 for good pixel behavior estimation. All the parameter values are empirically determined, and those with units are highly influenced by the calibration phantom material, water. Data processing and image reconstruction '. The experimental setup replaces all zerocounts in the measured data to 0.3 to reduce zero-count bias in line integrals (i.e. 6=0.3x1200 in Eq. (1), where 1200 is the projections per second.). Zero-count correction is described in D. Bushe, R. Zhang, G. H. Chen, and K. Li, “Unbiased zero-count correction method in low-dose high-resolution photon counting detector CT,” Phys. Med. Biol., vol. 68, no. 11, 2023, incorporated herein by reference in its entirety. Also, all measured counts (#) are converted to count-rates (# / s) by dividing them by the acquisition time for each projection. For image reconstruction, data from each detector row is reconstructed separately using a fan-beam reconstruction algorithm with a Hamming filter. The array size of the reconstructed slice is 950x950x1, and the voxel size is [0.5, 0.5, 0.2] mm. Subsequently, 24 adjacent slices are averaged to make an image with a slice thickness of ~5 mm.

[0233] Calibration methods'. The performance of three calibration methods are compared with the uncorrected results, as follows: (a) the multi-point flat field correction, which is called “MPFFC”; (b) the pixel-to-pixel variation correction using Eq. (14), called “PPVC”; and (c) the water linearization using Eq. (16), called “WL”.

[0234] (2) Assessment

[0235] In the experimental setup, ring and cupping / doming artifacts are assessed in reconstructed images of the test phantoms. Visual inspection, radial mean profile analysis, and radial mean standard deviation (RMSD) analysis, provide a comprehensive evaluation of both visual and quantitative aspects. Visual inspection offers an overall assessment, while the radial mean profile facilitates easier comparison between results. For the radial mean profile, first transform a reconstructed image from Cartesian to polar coordinates and then compute the mean values along the angular axis to reduce noise. Lastly, RMSDs provide a quantitative comparison. Two types of RMSDs are calculated to quantify ring artifacts (high frequency) and cupping / doming artifacts (low frequency). For the high-frequency RMSD, low-frequency components are removed by applying a high-pass filter to the image in polar coordinates. The high-frequency RMSD is defined as follows:

[0236] RMSDH= I^S^(gS(p) - aS)2(17)

[0237]

[0238]

[0239] where - (S * w)(p,< / >)} is the high-pass filtered image

[0240] averaged along the angular axis in polar coordinates and gp=p^ _pJjw(p). Pe{1< 2,

[0241] Pmax} is the radial index and < > G {1, 2, <pmax}

[0242] is the angular index. (p1and (p2represent the minimum and maximum angular indices used i for averaging. Here, g p, <p) is the original image in polar coordinates,

[0243]

[0244] and w(p, < / >) = (-P 2+,a2 )

[0245] e 2<r2is the Gaussian kernel. p1and p2represent the minimum and maximum radial indices used for RMSD calculation. In this study, we set the value of <7 to 50 and used a kernel window size of 101. The low-frequency RMSD is calculated as follows:

[0246] RMSDL(fi1,p2) = (18)

[0247]

[0248] y P2~P1! 1

[0249] where g

[0250]

[0251] ^(p) = * wKp>4>) is the average of the low pass filtered image in

[0252] polar coordinates and g^ =p> ^p2=Plg™(p)- For the low-frequency RMSD, the same

[0253]

[0254] Gaussian kernel is used as for the high-frequency RMSD.

[0255] FIG. 7 illustrates an example of a photon-counting Computed Tomography (PCCT) scanner system.

[0256] The calibration approach can be implemented in a photon-counting CT scanning system as described below with reference to FIG. 7. The X-ray CT apparatus 1 shown in FIG. 7 includes a gantry 10, a bed 30, and a console 40 that implements the processing of the medical imaging processing apparatus. In the present embodiment, the rotation axis of a rotation frame 13 in the non-tilted state, or the longitudinal direction of a tabletop 33 of the bed 30, is defined as a “Z-axis direction;” the axial direction orthogonal to the Z-axis direction and horizontal to the floor is defined as an “X-axis direction;” and the axial direction orthogonal to the Z-axis direction and vertical to the floor is defined as a “Y-axis direction.”

[0257] In an embodiment, the gantry 10 and the bed 30 are installed in a CT examination room, and the console 40 is installed in a control room adjacent to the CT examination room. The console 40 is not necessarily installed in the control room. For example, the console 40 can be installed together with the gantry 10 and the bed 30 in the same room. In any case, the gantry 10, the bed 30, and the console 40 are communicably connected to one another by wire or radio.

[0258] The gantry 10 is a scanner with a configuration for performing X-ray CT imaging on a subject (or an imaging object) P. The gantry 10 includes an X-ray tube 11, an X-ray detector 12, a rotation frame 13, an X-ray high voltage device 14, a controller 15, a wedge filter 16, a collimator 17, and a data acquisition system (DAS) 18.

[0259] The X-ray tube 11 is a vacuum tube that generates X-rays by emitting thermal electrons from the cathode (filament) to the anode (target) in response to application of a high voltage and supply of a filament current from the X-ray high voltage device 14. Specifically, X-rays are generated by the thermal electrons colliding with the target. Examples of the X-ray tube 11 include a rotating anode type X-ray tube that generates X-rays by emitting thermal electrons to the rotating anode. The X-rays generated in the X-ray tube 11 are, for example, formed into a cone-beam shape by the collimator 17, and applied to the subject P.

[0260] The X-ray detector 12 detects X-rays that have been emitted by the X-ray tube 11 and have passed through the subject P, and outputs an electrical signal corresponding to the X-ray dose to the DAS 18. The X-ray detector 12 includes a plurality of X-ray detection element lines, each including a plurality of X-ray detection elements aligned in the channel direction (the X-axis direction, or the column direction) along an arc having a center at the focus of the X-ray tube 11, for example. The X-ray detector 12 has an array structure in which a plurality of X-ray detection element lines, each including a plurality of X-ray detection elements aligned in the channel direction, are aligned in the segment direction (the Z-axis direction, or the row direction).

[0261] Specifically, the X-ray detector 12 can be, for example, a direct conversion type detector including a semiconductor element that converts incident X-rays into an electrical signal. The X-ray detector 12 is an example of the PCD according to the present embodiment, and will also be referred to as a “PCD 12.”

[0262] The rotation frame 13 supports an X-ray generator and the X-ray detector 12 rotatably around a rotation axis. Specifically, the rotation frame 13 is an annular frame that supports the X-ray tube 11 and the X-ray detector 12 in such a manner that the X-ray tube 11 faces the X-ray detector 12, and rotates the X-ray tube 11 and the X-ray detector 12 under the control of a controller 15 to be described later. The rotation frame 13 is rotatably supported by a stationary frame (not shown) made of a metal such as aluminum. Specifically, the rotation frame 13 is connected to an edge portion of the stationary frame via a bearing. The rotation frame 13 rotates around the rotation axis Z at a predetermined angular velocity while receiving power from a driver of the controller 15.

[0263] In addition to the X-ray tube 11 and the X-ray detector 12, the rotation frame 13 includes and supports the X-ray high voltage device 14 and the DAS 18. Such a rotation frame 13 is housed in an approximately-cylindrical case with a bore 19 constituting an imaging space. The bore approximately corresponds to the FOV. The central axis of the bore corresponds to the rotation axis Z of the rotation frame 13. Detection data generated by the DAS 18 is transmitted, for example, from a transmitter (not shown) to a receiver (not shown) arranged on a non-rotating portion (such as the stationary frame, illustration omitted in FIG. 13) of the gantry, and then transferred to the console 40.

[0264] The X-ray high voltage device 14 includes: a high voltage generator including electrical circuitry such as a transformer, a rectifier, etc. and having the function of generating a high voltage to be applied to the X-ray tube 11 and a filament current to be supplied to the X-ray tube 11; and an X-ray controller configured to control an output voltage in accordance with the X-rays emitted by the X-ray tube 11. The high voltage generator can be of a transformer type, or an inverter type. The X-ray high voltage device 14 may be provided in the rotation frame 13 to be described later, or in the stationary frame (not shown) of the gantry 10.

[0265] The controller 15 includes processing circuitry including a central processing unit (CPU), etc., and a driver such as a motor or an actuator, etc. The processing circuitry includes, as hardware resources, a processor, such as a CPU or a micro processing unit (MPU), and a memory, such as a read only memory (ROM) or a random access memory (RAM). The controller 15 can be realized by an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), or another complex programmable logic device (CPLD) or simple programmable logic device (SPLD). The controller 15 controls the X-ray high voltage device 14 and the DAS 18, etc. in accordance with instructions from the console 40. The processor implements the above control by reading and executing a program stored in the memory.

[0266] The CPU can execute a computer program including a set of computer-readable instructions that perform the functions described herein, and the program is stored in any of the above-described non-transitory electronic memories and / or a hard disk drive, CD, DVD, FLASH drive or any other known storage media. Further, the computer-readable instructions may be provided as a utility application, background daemon, or component of an operating system, or combination thereof, executing in conjunction with a processor and an operating system known to those skilled in the art. Further, the CPU can be implemented as multiple processors cooperatively working in parallel to perform the instructions.

[0267] The controller 15 also has the function of performing operation control of the gantry 10 and the bed 30 in response to an input signal from an input interface 43 to be described later attached to the console 40 or the gantry 10. For example, the controller 15 performs control to rotate the rotation frame 13, control to tilt the gantry 10, or control to operate the bed 30 and the tabletop 33 in response to an input signal. The control to tilt the gantry 10 is implemented by the controller 15 rotating the rotation frame 13 around an axis parallel to the X-axis direction, based on tilt angle information input through the input interface 43 attached to the gantry 10. The controller 15 may be provided either in the gantry 10 or in the console 40. The controller 15 may be configured by directly integrating a program in the circuitry of the processor, instead of storing a program in the memory. In this case, the processor implements the above-described control by reading and executing the program integrated in the circuitry.

[0268] The wedge filter 16 is a filter for adjusting the dose of X-rays emitted from the X-ray tube 11. Specifically, the wedge filter 16 is a filter that allows X-rays emitted from the X-ray tube 11 to pass therethrough, and attenuates the X-rays so that the X-rays emitted from the X-ray tube 11 to the subject P exhibit predetermined distribution. For example, the wedge filter 16 (or bow-tie filter) is a filter obtained by processing aluminum so that it has a predetermined target angle and a predetermined thickness.

[0269] The collimator 17 is lead plates or the like for narrowing the application range of X-rays that have passed through the wedge filter 16, and includes a slit formed by combining the lead plates or the like. The collimator 17 may be referred to as an “X-ray diaphragm.” The DAS 18 generates digital data indicating counts of X-rays detected by the X-ray detector 12 (also referred to as “detection data”) for each of a plurality of energy bands (referred to as “energy bins” or simply as “bins”). The detection data is a set of a channel number and row number of a source X-ray detection element, a view number indicating a collected view (also referred to as a projection angle), and data of the count value identified by the energy bin number. The DAS 18 is implemented by, for example, an application specific integrated circuit (ASIC) on which a circuit element capable of generating detection data is mounted. The detection data is transferred to the console 40.

[0270] The bed 30 is a device to place thereon the subject P to be scanned and move the subject P, and includes a base 31, a bed actuator 32, a tabletop 33, and a support frame 34.

[0271] The base 31 is a case that supports the support frame 34 movably in the vertical direction.

[0272] The bed actuator 32 is a motor or actuator that moves the tabletop 33 on which the subject P is placed in the longitudinal direction of the tabletop 33. The bed actuator 32 moves the tabletop 33 in accordance with control by the console 40 or control by the controller 15. For example, the bed actuator 32 moves the table top 33 in the direction orthogonal to the subject P so that the body axis of the subject P placed on the table top 33 matches the central axis of the bore of the rotation frame 13. The bed actuator 32 may also move the tabletop 33 in the body axis direction of the subject P in accordance with X-ray CT imaging performed using the gantry 10. The bed actuator 32 generates power by driving at a rotation speed corresponding to the duty ratio of the drive signal from the controller 15. The bed actuator 32 is implemented by a motor, such as a direct drive motor or a servo motor.

[0273] The tabletop 33 provided on the top surface of the support frame 34 is a plate on which the subject P is placed. The bed actuator 32 may move not only the tabletop 33 but the support frame 34 in the longitudinal direction of the tabletop 33. The console 40 includes a memory 41, a display 42, an input interface 43, and processing circuitry 44. Data communication between the memory 41, the display 42, the input interface 43, and the processing circuitry 44 is performed via a bus. The console 40 is described as being separate from the gantry 10, but the gantry 10 may include the console 40 or part of each constituent element of the console 40.

[0274] The memory 41 is a storage device, such as a hard disk drive (HDD), a solid state drive (SSD), or an integrated circuit storage device, etc., which stores various types of information. The memory 41 stores, for example, projection data and reconstructed image data. The memory 41 may be not only the HDD, SSD, or the like, but a driver that writes and reads various types of information in and from, for example, a portable storage medium such as CD, DVD, or a flash memory, or a semiconductor memory such as a random access memory (RAM). The storage area of the memory 41 may be in the X-ray CT apparatus 1, or in an external storage device connected via the network. For example, the memory 41 stores data of a CT image or a display image. The memory 41 also stores a control program according to the present embodiment.

[0275] The display 42 displays various types of information. For example, the display 42 outputs a graphical user interface (GUI) or the like for receiving a medical image (CT image) generated by the processing circuitry 44, and various types of operations from the operator. For the display 42, for example, a liquid crystal display (LCD), a cathode ray tube (CRT) display, an organic electro luminescence display (OELD), a plasma display, or any other display can be used as appropriate. The display 42 may be provided in the gantry 10. The display 42 may either be a desktop type or configured by a tablet device capable of wirelessly communicating with the console 40.

[0276] The input interface 43 receives various types of input operations from the operator, converts a received input operation into an electrical signal, and outputs the electrical signal to the processing circuitry 44. For example, the input interface 43 receives, from the operator, a collection condition for collecting projection data, a reconstruction condition for reconstructing a CT image, and an image-processing condition for generating a postprocessing image from the CT image, etc. For the input interface 43, for example, a mouse, a keyboard, a trackball, a switch, a button, a joystick, a touch pad, or a touch panel display can be used as appropriate. In the present embodiment, the input interface 43 does not necessarily include a physical operation component such as a mouse, a keyboard, a trackball, a switch, a button, a joystick, a touch pad, or a touch panel display. For example, the input interface 43 also includes electrical signal processing circuitry that receives an electrical signal corresponding to an input operation from an external input device provided separately from the apparatus, and outputs the electrical signal to the processing circuitry 44. The input interface 43 may be provided in the gantry 10. The input interface 43 may be configured by a tablet device capable of wirelessly communicating with the console 40.

[0277] The processing circuitry 44 controls the overall operation of the X-ray CT apparatus 1 in accordance with the electrical signal of the input operation output from the input interface 43. For example, the processing circuitry 44 includes, as hardware resources, a processor such as a CPU, an MPU, or a graphics processing unit (GPU), and a memory such as a ROM or a RAM. With a processor that executes a program loaded into the memory, the processing circuitry 44 performs a system control function 441, a pre-processing function 442, a reconstruction function 443, and a display control function 444. Each of the functions (the system control function 441, the pre-processing function 442, the reconstruction function 443, and the display control function 444) is not necessarily implemented by a single processing circuit. Processing circuitry can be configured by combining a plurality of independent processors, and the processors can execute respective programs to implement the functions. The system control function 441 controls each function of the processing circuitry 44 based on an input operation received from the operator via the input interface

[0278] 43. Specifically, the system control function 441 reads a control program stored in the memory 41, loads it into a memory in the processing circuitry 44, and controls each part of the X-ray CT apparatus 1 in accordance with the loaded control program. For example, the processing circuitry 44 performs each function of the processing circuitry 44 based on an input operation received from the operator via the input interface 43. For example, the system control function 441 obtains a two-dimensional positioning image of the subject P to determine the scan range, imaging condition, etc. The positioning image can also be referred to as a “scanogram” or “scout image.”

[0279] The pre-processing function 442 generates data obtained by performing preprocessing on detection data output from the DAS 18. Data (detection data) before preprocessing and data after pre-processing can be collectively referred to as “projection data.” The pre-processing function 442 is an example of the pre-processor.

[0280] The reconstruction function 443 generates CT image data by performing reconstruction processing using a filtered back projection method, a successive approximation reconstruction method, a stochastic image reconstruction method, or the like, on the projection data generated by the pre-processing function 442. The reconstruction function 443 is an example of the reconstruction processor. Image filtering, smoothing, volume rendering, or image differential processing can be applied to the CT image data if required. The display control function 444 converts CT image data generated by the reconstruction function 443 into tomographic image data of a given cross section, or three-dimensional image data by a publicly-known method, based on the input operation received from the operator via the input interface 43. The generation of three-dimensional image data can be performed directly by the reconstruction function 443. The display control function 444 is an example of the display controller.

[0281] In one implementation, the X-ray tube 11 is a single source emitting a broad spectrum of X-ray energies, and the PCD 12 can use a direct X-ray radiation detectors based on semiconductors, such as cadmium telluride (CdTe), cadmium zinc telluride (CZT), silicon (Si), mercuric iodide (HgI2), and gallium arsenide (GaAs). As mentioned above, semiconductor-based direct X-ray detectors generally have much faster time response than indirect detectors, such as scintillator detectors. The fast time response of direct detectors enables them to resolve individual X-ray detection events, although at the high X-ray fluxes typical in clinical X-ray applications, some pileup of detection events may occur. The energy of a detected X-ray is proportional to the signal generated by the direct detector, and the detection events can be organized into energy bins yielding spectrally resolved X-ray data for spectral CT.

[0282] RESULTS

[0283] (1) Calibration results

[0284] FIG. 8A shows the optimal location parameters from Step 1 of the calibration process. The different shapes represent different input offsets. Phantom sizes include 18 cm, 24 cm,

[0285] (= / I(rx)2+ ( (ry ) )wr’teninpl°taresimilar to their

[0286]

[0287] input offsets (digital positioning values). The goodness of the location parameter estimation can be visually assessed by collecting data from different CT scans into a single graph. In FIG. 8B, data collected from different offsets and phantom sizes are plotted for a central pixel. The alignment between data clusters is crucial, as deviations could cause calibration inaccuracies in the following calibration steps. The data clusters from the estimated location parameters show close alignment. FIG. 9 shows the curve fitting results of count-rates against water thickness at tube current 200 mA, derived from the second step of the calibration process. Results are presented for total count, bin 1, and bin 5 of the peripheral pixel and the center pixel. Due to the lower beam intensity in the peripheral pixels, their data are noisier than those of the central pixels, resulting in a greater spread along the vertical axis. The data discontinuity observed in the thick water region of the peripheral pixel plots is because counts in those areas are close to zero. The total count plots show relatively straight fitting lines for both the central pixel and the peripheral pixel, while the bin count plots exhibit slight curvatures. Notably, the bin 1 and bin 5 plots of the center pixel show different curvature at high log count-rates, likely caused by the pulse pileup effect. Despite the significantly different conditions of the data given, the proposed curve-fitting method demonstrates good fitting for different count data in both pixels.

[0288] FIG. 10 shows the good pixel behavior fitting results for water thicknesses of 0 cm and 20 cm, obtained from the third step of the calibration process. The count-rates across pixels exhibit bell-shaped distributions, which is largely due to the bowtie filter. For both total count and different energy bins, the fitting results demonstrate smooth bell curves across all pixels and properly passes through the center region of pixel variations.

[0289] (2) Correction results

[0290] FIG. 11 presents reconstructed total count and energy bin images (bins 1-5) of 24 cm water phantoms without offset, generated using different calibration methods. The uncorrected images consistently show severe ring artifacts, but different ring artifact patterns are observed across different energy bins. Cupping and doming artifacts are also observed, with cupping and doming appearing in the center of the bins 1-3 images and the bin 5 image, respectively. The PPVC images show outstanding ring artifact reduction for both total count and all energy bins. Cupping and doming artifacts in the center of the bin images, however, seemingly remain similar to the uncorrected images, appearing cupping in bins 1- 3 and doming in bin 5. The WL images also show outstanding performance in ring artifact reduction across both total count and all energy bins, with significant improvement in the cupping and doming artifacts compared with those in the PPVC images. The MPFFC images overall show a decent ring artifact reduction for both total count and all energy bins, but they introduce more prominent doming artifacts for bins 1-3 and cupping artifacts for bin 5 than those in the uncorrected images. These findings are representative of trends observed in subsequent results.

[0291] FIG. 12 presents reconstructed total count and energy bin images of 24 cm water phantoms with a 5 cm offset. Overall, all the images exhibit similar levels of ring artifacts as in FIG. 11. While the uncorrected and PPVC images still show these cupping / doming artifacts in bins 1-3 and bin 5, they are now centered around the FOV center rather than the phantom's center. In contrast, the MPFFC images maintain cupping and doming artifacts centered around the phantom's center.

[0292] In FIG. 13, radial mean profiles are presented for the total count, bin 1, and bin 5 images in FIGs. 11 and 12. The radial mean profiles were calculated based on the circular sectors in each reconstructed image in FIGs. 11(1 a) and 12(1 a), respectively. The circular sectors indicate the angular averaging ranges, defined by the span between central angles

[0293]

[0294] and 2- All the profiles of the uncorrected images 1302 show high levels of rings. The MPFFC profiles 1308 show reduced rings compared to the uncorrected profiles, but we can observe more noticeable cupping / doming artifacts. The profiles of the PPVC 1304 and WL 1306 images show remarkably reduced rings compared to the uncorrected and the MPFFC profiles. Notably, the total count profiles of the PPVC and WL are nearly identical and consistent across the entire radial range for both the no-offset and 5 cm offset results. In all the plots, the PPVC profiles 1304 consistently pass through the center of the uncorrected profiles. The WL profiles 1306 generally show better flatness compared to the PPVC profiles, although some residual cupping and doming effects remain in the bin 1 and bin 5 profiles.

[0295] FIG. 14 presents high- and low-frequency RMSDs of the total count and bin 1-5 images in FIGs. 11 and 12 for the different two offsets. The RMSDs were calculated based on the yellow circular sectors in the reconstructed images in FIG. 1 l(la) and FIG. 12(la). Both the high-frequency and low-frequency RMSDs of each method present very similar results across the two offsets. The MPFFC 1408 presents high-frequency RMSDs that are over three times smaller than those of the uncorrected 1402 results across all cases. Both the PPVC 1404 and the WL 1406 show significantly greater improvement than the MPFFC 1408 for all cases and are closely aligned with each other. For total counts, the MPFFC 1408, PPVC 1404, and WL 1406 improve high-frequency RMSDs by 88.4%, 95.4%, and 95.7%, respectively, for no offset, and by 87.2%, 95.2%, and 95.5%, respectively, for 5 cm offset. For energy bin counts, they achieve average improvements across all energy bins of 74.5%, 92.5%, and 93.6%, respectively, for no offset, and by 76.2%, 92.4%, and 93.5%, respectively, for 5 cm offset.

[0296] For low-frequency RMSDs, the MPFFC 1408 show even worse RMSDs than the uncorrected results for bins 1-2 at both offsets and bin 3 for no offset. Generally, the PPVC 1404 shows improved RMSDs over the uncorrected results. The WL 1406 generally shows improvements over the PPVC 1404, especially for bin 5. The MPFFC 1408 shows average enhancements of 65.3% for total counts and -95.1% for energy bin counts. For total counts, the MPFFC 1408, PPVC 1404, and WL 1406 improve low-frequency RMSDs by 49.8%, 93.3%, and 97.2%, respectively, for no offset, and by 50.6%, 97.4%, and 97.9%, respectively, for 5 cm offset. For energy bin counts, they achieve average improvements across all energy bins of -120.0%, 54.4%, and 55.7%, respectively, for no offset, and by -27.1%, 47.3%, and 70.4%, respectively, for 5 cm offset.

[0297] FIG. 15 shows reconstructed head phantom images, with a narrower display window applied to the head images: [-100 100] HU. The PPVC and WL effectively reduce ring artifacts, with WL also significantly minimizing cupping and doming artifacts. The MPFFC reduces ring artifacts, especially for total counts, but artifacts persist in energy bin images along with cupping / doming effects.

[0298] FIG. 16 shows reconstructed body phantom images. Both the PPVC and WL show remarkable improvements in ring artifact reduction in the primary phantom material. Also, it is important to note that they effectively reduce ring artifacts of different materials in the concentration inserts. Unlike other phantom results, the WL introduces slight doming artifacts across all energy bins, but not in the total count. The MPFFC effectively improves ring artifacts but still introduces severe doming artifacts in bins 1-3 and cupping artifacts in bin 5. Additionally, no method effectively corrects beam-hardening artifacts, around the calcium 300 mg / ml insert (see arrows in FIG. 16(3c)), particularly in the lower energy bin images.

[0299] A PCD calibration method is disclosed that addresses ring artifacts in energy bin images caused by pixel-to-pixel variations, as well as count-rate non-linearity related to water thickness. Findings demonstrate that the calibration process significantly enhances image quality across a variety of phantoms. The WL improved high-frequency RMSDs by 95.5-95.7% (total counts) and 93.5-93.6% (energy bins), and low-frequency RMSDs by 97.2-97.9% (total counts) and 55.7-70.4% (energy bins), showing robust artifact reduction.

[0300] Doming artifacts are observed across all energy bins in the body phantom WL results. It is mainly due to zero-count bias, which persists despite the zero-replacement (zeroreplacement with 0.3 introduce a positive post-log bias.). FIG. 17A shows the average difference of ROI 2 relative to ROI 1 in the WL results (see FIG. 6B for ROI locations). The ROI 2 averages for all energy bins are consistently lower than those of ROI 1. In contrast, total count and total count with bin correction averages are consistent between the ROIs, indicating that individual bin corrections match total count corrections, with discrepancies arising within individual bins. FIG. 17B highlights zero-count occurrence as red points in sinograms. Zero-count occurrence rates in the body phantom were 0.35-1.42% for bins 1-5, and 0.02% for total count.

[0301] The results reveal bin-dependent issues that vary across the calibration methods. The MPFFC generally introduces more pronounced cupping / doming artifacts compared to the uncorrected images, with these artifacts centered on the phantom rather than the FOV center. This is likely due to spectral discrepancies between the calibration and test data. As X-ray spectra are hardened by attenuation through the imaged object, the differences between the air scan and object scan spectra become more pronounced, causing air scan-based coefficients to introduce unintentional weightings when applied to object-attenuated data. These weightings can either amplify or reduce existing low-frequency artifacts, depending on the energy bin, resulting in distinct low-frequency RMSDs that differ from those produced by water-based calibration methods. In contrast, the PPVC maintained similar levels of cupping and doming artifacts as seen in the uncorrected images, demonstrated by the radial mean profiles in Fig. 13, confirming effective good-pixel correction. The WL method provided relatively consistent correction effects for these artifacts across energy bins, showing the lowest bin-dependency.

[0302] The MPFFC shows relatively higher high-frequency RMSDs at bin 5, contrasting with the results from the proposed calibration methods. It is likely due to increased count-rate non-linearity against water thickness caused by the pulse pileup effect. For calibration methods that use multiple data points, like the MPFFC, capturing higher degrees of non- linearity is challenging unless additional data conditions are acquired, as these methods approximate coefficients for intermediate conditions between existing calibration conditions. In contrast, the proposed calibration methods show relatively consistent high-frequency RMSDs across all energy bins, proving their capability to account for higher degrees of nonlinearity. This is achieved by acquiring near-continuous thickness of the calibration phantom through CT scans.

[0303] The disclosed method uses only a single material, water, for calibration, which may lead to the remaining beam-hardening artifacts caused by materials with higher atomic numbers, as observed in the WL images in FIG. 16. Increasing the number of calibration materials, however, makes it impossible to handle each energy bin independently without additional image processing, such as segmentation, necessitating MD that utilizes multiple energy bins simultaneously. While these integrated calibration methods provide distinct advantages, the disclosed method offers flexibility by enabling the resolution of pixel-to-pixel variation as an independent issue.

[0304] This flexibility allows the disclosed method to support several potential approaches to address the remaining beam-hardening artifacts: MD following PPVC, MD following WL, or two-pass BHC (high-Z BHC following WL). Using either PPVC or WL, the proposed method enables subsequent MD through a simple calibration process, assuming uniformity across all pixels. This MD can be performed in either the projection or the image domain. The high-Z BHC in the two-pass BHC approach specifically targets high-Z material beamhardening effects without requiring multiple energy bins. Instead, it includes an additional segmentation step in the image domain to estimate the thickness and density of high-Z materials. See P. M. Joseph and R. D. Spital, “A Method for Correcting Bone Induced Artifacts in Computed Tomography Scanners,” J. Comput. Assist. Tomogr., vol. 2, no. 1, pp.

[0305] 100-108, 1978, incorporated herein in its entirety. In the disclosed calibration method, phenomena such as charge sharing, K-fluorescence, and pulse pileup were not explicitly modeled. These effects may be partially addressed through the calibration process, as it is based on comprehensive measurements that inherently contain combined effects; however, they are not entirely resolved. The absence of models for these factors introduces limitations to the proposed method’s performance, especially in cases where the test data exhibit effects that significantly differ from those in the calibration data. Incorporating models for these effects into the calibration process could enhance the robustness of the calibration method by providing more accurate corrections. The proposed method, however, emphasizes prioritizes simplicity and ease of implementation, serving as a foundational framework for further development. Future studies can build on this approach by integrating models for charge sharing, K-fluorescence, and pulse pileup, further improving its accuracy and broadening applicability across diverse imaging tasks.

[0306] The disclosure presents a calibration framework tailored for PCDs to address pixel-to-pixel variations and minimize ring artifacts in CT images. The disclosed framework incorporates a streamlined calibration data acquisition strategy using CT scans of calibration phantoms, coupled with robust calibration steps. The disclosed method demonstrated significant improvements in both high-frequency and low-frequency RMSD reduction, effectively addressing ring and cupping / doming artifacts, respectively. It outperformed a conventional multi-point calibration method, especially in handling the complex non-linear behavior of PCDs. Overall, the proposed calibration framework provides a robust and flexible solution for optimizing PCD performance in CT imaging.

Claims

CLAIMS:

1. A method for calibrating a photon-counting computed tomography (PCCT) apparatus, the method comprising:performing, using the PCCT apparatus, for each offset position of a plurality of offset positions from a rotation center, a scan of a calibration phantom placed at the offset position to obtain calibration data including measured count rate data, wherein the PCCT apparatus is configured with a semiconductor-based photon-counting detector (PCD) having a plurality of pixels;estimating, by processing circuitry, for each of the plurality of offset positions, location parameters of the calibration phantom based on the obtained calibration data;determining, by the processing circuitry based on the estimated location parameters, thickness data of the calibration phantom;estimating, by curve fitting using the measured count rate data and the determined thickness data, a curve indicating a relationship between the energy bin count rate and phantom thickness, for each energy bin of each pixel of the plurality of pixels;generating, by the processing circuitry based on the estimated curves of each energy bin of each pixel, a first look-up table for pixel-to-pixel count rate variation correction; and applying, by the processing circuitry, the first look-up table to acquired energy bin count data to minimize artifacts in an image reconstructed using the acquired energy bin count data.

2. The method of claim 1, wherein the step of estimating the location parameters further comprises creating a three-dimensional image of the calibration phantom, as a template image, using information about a geometry of the PCCT apparatus and a shape of the calibration phantom.

3. The method of claim 2, wherein the step of estimating the location parameters further comprises creating, by the processing circuitry, a weight image that assigns higher importance to certain regions in the template image.

4. The method of claim 3, wherein the step of estimating the location parameters further comprises optimizing, by the processing circuitry, an alignment between the template image and the calibration data by adjusting a position of the template image.

5. The method of claim 4, further comprising iteratively estimating, by the processing circuitry, the plurality of offsets of the template image until the template image substantially matches the calibration data so as to determine the location parameters.

6. The method of claim 5, further comprising:determining, by the processing circuitry, a quality of the alignment using a normalized cross-correlation (NCC) process, wherein the NCC process determines how similar the template image is to the calibration data, ignoring overall brightness differences; and after the optimal alignment is determined, estimating, by the processing circuitry, the thickness data of the calibration phantom based on the location parameters.

7. The method of claim 1, wherein the step of estimating the curve further comprises: applying a curve-fitting technique to obtain a smooth curve through data points of the measured count rate data points;estimating a count rate at spaced thicknesses using interpolation; andwhen a thickness point exceeds a measured range, extending the curve to determine the count rates for those thickness points, using extrapolation.

8. The method of claim 1, wherein the step of generating the first look-up table comprises:estimating a behavior at an ideal pixel according to a predetermined model; calculating pixel correction coefficients based on expected behavior of the ideal pixel so that each pixel matches the behavior of the ideal pixel; andstoring the calculated pixel correction coefficients in the first look-up table.

9. The method of claim 1, further comprising:generating, by the processing circuitry, a second look-up table for pixel-to-pixel count rate variation correction with respect to measurement of water; andapplying, by the processing circuitry, the first and second look-up tables to the acquired energy bin count data to minimize the artifacts.

10. The method of claim 9, wherein the step of generating the second look-up table further comprises:determining water linearity correction coefficients such that all of the pixels behave in a same way when measuring water; andstoring the determined water linearity correction coefficients in the second look-up table.

11. The method of claim 1, further comprising correcting, by the processing circuitry using the first and the second look-up tables, PCCT scan data to obtain a corrected image, in which all parts of the corrected image have a uniform response.

12. The method of claim 11, wherein the applying step further comprises determining, using the first lookup table, a pixel correction coefficient for an actual count rate of the acquired energy bin count data using linear interpolation.

13. The method of claim 12, wherein the applying step further comprises calculating, by the processing circuitry, an interpolation ratio given a measured line integral for the acquired energy bin count data.

14. The method of claim 12, wherein the applying step further comprises performing, by the processing circuitry, pixel-to-pixel variation correction and water linearization correction sequentially.

15. The method of claim 13, further comprising calculating, by the processing circuitry, the interpolation ratio based on the measured line integral.

16. The method of claim 13, further comprising calculating, by the processing circuitry, the interpolation ratio based on the measured count-rate.

17. The method of claim 12, further comprising performing, by the processing circuitry, direct pixel-to-pixel variation conversion to line integrals.

18. The method of claim 9, further comprising:generating, by the processing circuitry, corresponding first and second look-up tables for discrete mA values and respective specific count-rate intervals for performing tube current modulation;for given test values, and for each look-up table of the first and the second look-up tables,identifying, by the processing circuitry, nearest neighboring mAs and countrates in each look-up table; andcomputing, by the processing circuitry, correction coefficients by combining four neighboring coefficients corresponding to the given test values.

19. An apparatus for calibrating a photon-counting computed tomography (PCCT) device, the apparatus comprising:processing circuitry configured tocause the PCCT apparatus to perform, for each offset position of a plurality of offset positions from a rotation center, a scan of at least one calibration phantom placed at the offset positions to obtain calibration data including measured count rate data, wherein the PCCT apparatus is configured with a semiconductor-based photon-counting detector (PCD) having a plurality of pixels,estimate, for each of the plurality of offset positions, location parameters of the at least one calibration phantom based on the obtained calibration data;determine based on the estimated location parameters, thickness data of the at least one calibration phantom;estimate, by curve fitting using the measured count rate data and the determined thickness data, a curve indicating a relationship between the energy bin count rate and phantom thickness, for each energy bin of each pixel of the plurality of pixels;generate, based on the estimated curves of each energy bin of each pixel, a first look-up table for pixel-to-pixel count rate variation correction;apply the first look-up table to acquired energy bin count data to minimize artifacts in an image reconstructed using the acquired energy bin count data.

20. The apparatus of claim 19, wherein the processing circuitry is further configured to:cause the PCCT apparatus to perform, for each offset position of the plurality of offset positions from a rotation center, a scan of a plurality of calibration phantoms, of multiple sizes, placed at the offset position to obtain calibration data including the measured count rate data..

Images