- Email: [email protected]
Determination of a neutron beam divergence after the rocking curve concept using Richardson–Lucy’s unfolding algorithm
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Contents lists available at ScienceDirect
Nuclear Inst. and Methods in Physics Research, A journal homepage: www.elsevier.com/locate/nima
Determination of a neutron beam divergence after the rocking curve concept using Richardson–Lucy’s unfolding algorithm E.S. Souza a ,∗, M.I. Silvani b , G.L. Almeida b , R.T. Lopes a a
Universidade Federal do Rio de Janeiro, COPPE, Centro de Tecnologia, Cidade Universitaria Bloco G, Ilha do Fundao, 21945-970 Rio de Janeiro - RJ, Brazil Instituto de Engenharia Nuclear, Reator Argonauta - CNEN, Rua Helio de Almeida 75, Cidade Universitária, Ilha do Fundao Caixa Postal 68550 CEP 21941-972 Rio de Janeiro - RJ, Brazil b
ARTICLE
INFO
Keywords: Neutron beam divergence Rocking curve L/D Richardson–Lucy
ABSTRACT Thermal neutron radiographs acquired under high beam divergences suffer the impact of intense penumbrae degrading their final quality. As a divergence reduction is not always feasible, one possible alternative is an a posteriori treatment to restore the degraded images, such as the Richardson–Lucy — RL unfolding algorithm. Such a procedure requires the characterization of the spoiling agent Point Spread function – PSF in order to apply it in the inverse way. It can be deduced from the blur in the radiograph of a shielding blade edge. This blur depends upon the beam divergence and the object–detector gap. Due to the complex scattering of neutrons along a reactor channel, it is usual to express this divergence as the inverse of an L/D ratio. A novel approach based on the concept of Rocking Curve – RC, a term borrowed from the X-ray diffraction field, has been recently proposed which yielded slightly better quantitative results. After this concept, every sub-element of a surface source emits neutrons anisotropically following a bell-shaped profile. The RC angular semi-width incorporating neutron scattering and geometrical blur, is assigned as the beam divergence. The present work aims at its assessment through a quantitative determination ratified by a qualitative evaluation of radiographs unfolded by the RL algorithm. Its main purpose is an additional ratification of the soundness, consistency and robustness of the RC concept by comparison with formerly obtained quantitative results. In spite of the utterly different approaches and techniques, the outcome has corroborated the novel concept. All data treatment is simple and performed by an ad hoc written Fortran 90 program embedding the required algorithms.
1. Introduction
depends upon the beam divergence w, being amplified by the objectdetector gap g.
Radiographic images acquired with thermal neutrons are infested by plagues such as beam divergence, poor detector resolution, neutron scattering, statistical fluctuation and electronic noise which degrade their final quality. All of them are impossible to eliminate or at least hard to mitigate at their origin due to their fundamental features or technological and engineering constraints. Yet, although not possible to restore the primordial image, it is feasible to improve them through a post-acquisition treatment such as an unfolding with the Richardson– Lucy — RL algorithm [1,2]. This treatment requires however the characterization of the spoiling agent Point Spread Function–PSF in order to apply it in the inverse way. The PSF is the response of an image acquisition system to a pointlike input. For isotropic systems it has an intensity distribution with a radial symmetry engulfing all the spoiling agents. Although hard to characterize, it can be fairly expressed by the rotation of a Gaussian around its epicenter, and its standard deviation deduced from the blur cast on a detector by the straight edge of shielding blade. This blur
As for the beam divergence itself, it is worthwhile to define its meaning within the scope of this work as follows. A flat surface source emitting only in the perpendicular direction, should not produce any penumbra at the detector, independently of the object-detector gap or its thickness. As early stated by Berger [3] – quoted in Domanus [4] – ‘‘. . . the present state of art is such that parallel beams are definitely preferred’’. This assertion was later on refuted by Barton [5] who concluded that a divergent collimator produces better images. This disagreement seems to be a matter of misinterpretation. Perhaps Berger had in mind an ideal flat source as above described and not a real one where each source sub-element would emit neutrons isotropically. Or perhaps he was actually referring to quasi-parallel beams, as those tied to high L/D ratios. Under such a circumstance, a divergent collimator would naturally produce a better image, as stated by Barton, because in this case, the ratio of the source–detector clearance L to source size D
∗ Corresponding author. E-mail addresses: [email protected] (E.S. Souza), [email protected] (M.I. Silvani), [email protected] (G.L. Almeida), [email protected] (R.T. Lopes).
https://doi.org/10.1016/j.nima.2018.08.055 Received 15 May 2018; Received in revised form 2 August 2018; Accepted 18 August 2018 Available online xxxx 0168-9002/© 2018 Published by Elsevier B.V.
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
would increase towards that exhibited by an ideal point source, where the L/D ratio would obviously raise to infinite. An ideal point-source as well, would not blur the image, and thus is treated in this work as nondivergent, in spite of its conical beam. In a real neutron acquisition system however – where a virtual flat source could be visualized somewhere in the reactor neutron channel – scattering causes some neutrons to hit the detector at different angles, as if they were delivered from several sources at different distances and intensities. The overall impact on the quality of final image is equivalent to that caused by a single source of unknown size D located at an as well unknown distance L from the detector. As comprehensively addressed by Domanus [4], other techniques e.g., [6–10] rather than simple length measurements should be carried out to measure the effective L/D, which is always lower than the purely geometric one and rules the quality of the final image. Hence, disregarding the specific geometric arrangement at any facility, the images would exhibit equivalent qualities if they were acquired under the same effective L/D. The approach condensing the spoiling agents into a single L/D parameter overcomes the complex task of dealing with them, mainly with neutron scattering. Recent works have been proposed [11,12] positing the concept of Rocking Curve – RC, a designation borrowed from the X-ray diffraction field. It assumes that the source emit neutrons anisotropically with intensities ruled by a 3D-Gaussian distribution, and the beam divergence expressed by its angular half width at half maximum – HWHM designated as w in this work. For the sake of completeness, some features of those works including drawings, sketches and equations are shared by the present work. Although both earlier works share the same concept, they employ utterly different approaches and techniques. While Ref. [11] employs a position sensitive detector, a slit collimator and data treatment involving unfolding and extrapolation, Ref. [12] requires the neutron radiograph of a shielding blade edge, and an ad hoc written algorithm. Despite these differences, the results agree within 7% as follows: Reference RC HWHM (min)
[11] 75.91 ± 1.52
the searched one tied to the RC HWHM 𝑤0 . Further details are presented in Sections 2.4 and 2.5. An illustrative concept of the RC is sketched in Fig. 1, where a high divergence corresponds to a broad Gaussian (large w) while a zero divergence would be expressed by an infinitely narrow one, emitting neutrons solely in the axial direction. 2. Materials and methods The direct experimental determination of the PSF is very hard due to the low counting statistics imposed by its small size. Instead, it is replaced by the measurement of the Line Spread function – LSF, the response of a system to a line which can be obtained by differentiation of the Edge Response Function – ERF. This function is provided by the radiograph of a shielding blade edge. A Gaussian is then fitted to the LSF distribution and its FWHM assigned as the PSF width, since LSF and PSF share the same FWHM when the distribution follows a Gaussian profile [14] as depicted in Fig. 2. To overcome the related radiological burden, reactor operation costs and a cumbersome work, the family of curves s(w, g ) is obtained by an ad hoc written algorithm which generates synthetic radiographs and carries out the required treatments. 2.1. Generation of the synthetic images To achieve this task a virtual 2D detector emulating a real one such as an imaging plate is positioned at the end of a neutron channel perpendicularly to its axis, as sketched in Fig. 3. Somewhere in the neutron channel, a flat neutron source filling its whole cross-section would emit neutrons which could or not be intercepted by a shielding blade placed at a chosen distance g from the detector. Since its edge is aligned with the vertical axis dividing the image into two regions, the left region would be shielded while its right companion would be hit by the neutrons coming from the source. This ideal situation would occur only if the beam divergence or g were zero. For all other cases the penumbra would invade the neighborhood of the border line blurring the otherwise sharp edge. The intensity of the neutrons hitting the detector is assumed to depend upon the angle 𝛷 of their paths to the normal to the source as illustrated by the pictured bell-shaped surface. Once a spatial detector resolution 𝛿 is defined, an MxN matrix is assigned to the image, so that 𝛿.M and 𝛿.N represent its width and height respectively. Each pixel intensity is obtained by the summation of all neutrons coming from the source subdivided like the detector into KxH square elements of size 𝛼. So, 𝛼.K and 𝛼.H represent the width and height of the neutron channel and its virtual source. For a real source-detector set, the pixel intensity would depend upon the source intensity, detection efficiency and exposure time, but since that – for tiff images – it should not surpass the limit of 65,535, after summation, the pixel intensities as defined by Eq. (1) are normalized to this limit.
[12] 81.20 ± 1.18
Besides this fair agreement, the single blade required by Ref. [12], instead of an elaborated test-object to measure the L/D, makes the method affordable and attractive due to its simplicity. Furthermore, its basic algorithm used to measure the RC has been also used to measure the L/D agreeing fairly with Ref. [13] – based on neutron flux measurements – and exhibiting a lower uncertainty than the standard method presented in Ref. [6] as follows: Reference L/D ratio L/D uncertainty %
[6] 12.5–16.7
[12] 20.14 ± 1.17 5.8
[13] 19.2 ± 0.71 3.7
The above data refer to the neutron port of the Argonauta research reactor at Instituto de Engenharia Nuclear – CNEN, Rio de Janeiro Brazil, where the measurements [11–13] have been done. In spite of the better uncertainty exhibited by Ref. [13], it is less employed due to its timeconsuming. Taking into account this performance, and considering that the RC concept has not been yet comprehensively and exhaustively crosschecked, it is an advisable policy to ratify its soundness, consistency and robustness by another approach. Within this frame, the present work addresses the determination of the beam divergence through unfolding of a neutron radiograph with the RL algorithm. This is done by an analysis of radiographs unfolded with different PSF widths s, which selects the better amidst them. Its related PSF width 𝑠𝑥 and the object-detector gap 𝑔0 used to get the original radiograph constitute the coordinates [g0 , 𝑠x ] from which the related beam divergence 𝑤0 will be inferred. To accomplish this task, a family of synthetic curves s(w, g ) is generated and the curve s(𝑤0 , g) hit by [g0 , 𝑠x ] is assigned as
𝑝(𝑖, 𝑗) =
𝐾 ∑
𝛺(𝑖, 𝑘).
𝑘=1
𝛷𝑤 = 𝑤.][𝑙𝑛4
𝐻 ∑
{ [ ]2 } exp −0.5 𝛷 (𝑖, 𝑗, 𝑘, ℎ) ∕𝛷𝑤
(1)
ℎ=1
(2)
where: p(𝑖, 𝑗) = Pixel intensity on the detector: 𝑖 = 1 to M, 𝑗 = 1 to 𝑁. 𝛺(𝑖, 𝑘) = Bump function: 0 if the neutron hits the shielding, 1 otherwise. 𝛷(𝑖, 𝑗, 𝑘, ℎ) = Angle between the normal to the source and the straight line connecting the points (𝑖, 𝑗) on the detector and (𝑘, ℎ) on the source. 𝛷𝑤 = Angular standard deviation of the generatrix Gaussian, for the 3D Rocking Curve. w = Angular Half-Width at Half Maximum − HWHM of the Gaussian. K = No. of elements along the source width. H = No. of elements along the source height. M = No. of pixels along the image width on the detector. N = No. of pixels along the image height on the detector. 173
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 1. Concept of Rocking Curve as an anisotropic emission characterizing the beam divergence defined by its half width at half maximum – HWHM w.
Fig. 2. Response of an image acquisition system to a point and a line. Since the Point Spread Function – PSF is hard to measure experimentally it is replaced by the Line Spread Function – LSF. For a Gaussian distribution, they share the same FWHM. Fig. 3. Generation of synthetic images. Each pixel on the detector – one is highlighted – collects and integrates the radiation coming from all source elements. Unless blocked by the shielding blade, its intensity depends upon the angle 𝛷 of their paths to the normal to the source as illustrated by the pictured bell-shaped surface.
It should be pointed out that 𝛺(𝑖, 𝑘) does not depend upon j and h, because the shielding edge has a vertical orientation. As sketched in Fig. 4, which is fairly self-explanatory, the bell-shaped RC – represented by the rotation of a Gaussian around its axis – rules the flux of neutron emitted by each source element. Unlike an isotropic emission, the flux intensity depends upon the angle 𝛷(𝑖, 𝑗, 𝑘, ℎ) between the normal to the source and the straight line connecting the pixel (𝑖, 𝑗) on the detector and the element (𝑘, ℎ) on the source. Source and detector share the same dimensions of the neutron channel, but different pixel sizes. Both can be freely chosen as input data, but in order to emulate the detector actually used to acquire the experimental radiographs, a 0.05 mm resolution has been assigned to the virtual one. As for the source, the elements can be set as small as wished, provided that the choice does not enlarge the processing time to unreasonable values. Other parameters which should be entered as input data for the 3D analytical geometry treatment are w, g, L, and D.
As it is assumed that the RC can be represented by a Gaussian function, some limit should be set to discard insignificant contributions arising from its tails. However to keep a conservative approach, all contributions coming from angles below F.w are taken into account. An F -value of 5.0 has been employed in this work. 2.2. Determination of the ERF, LSF and PSF width for the synthetic images An ERF is obtained by averaging the pixel intensities along a line parallel to the edge shadow, both for synthetic and experimental images. This averaging encompassing several pixels along a row aims at to smooth the unavoidable noise and statistical fluctuation, and is written 174
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 4. 3D Rocking Curve as an extended concept borrowed from the X-ray diffraction field using a Gaussian as generatrix. The flux intensities 𝑆𝑎 and 𝑆𝑏 are tied to the respective angles 𝛷𝑎 and 𝛷𝑏 .
is worth to review briefly the mechanism of this algorithm as depicted in Figs. 7 and 8 due to its paramount role in this work. 2.3. Brief overview of the Richardson–Lucy unfolding algorithm In order to simplify the description it is convenient to choose a 1D example, such as a gamma-ray spectrum. In the absence of noise of any kind, and with an ideal zero-resolution detector, the gamma-rays would appear as vertical lines along the energy axis. If these conditions are not fulfilled, the contents (counts/channel) of every channel leaks to neighbor companions contaminating them with contribution of alien energies. In consequence of this convolution, the otherwise vertical lines are degraded to the usual peak shape. Yet, the information is not lost – just concealed – being capable to be retrieved by a proper treatment as schematized in Fig. 7. The RL algorithm employs an iterative procedure to determine the 𝑡 vector 𝑢𝑡+1 𝑗 from its predecessor 𝑢𝑗 through the correcting vector 𝑓𝑗 . Hence, when 𝑓𝑗 approaches the unitary vector, the spectra become more similar to each other after an asymptotical behavior. Expressing in words, what is already precisely described by the equations, the observed spectrum 𝑐𝑖 (right) is considered as if it were a non-degraded one. The vertical strips represent the primordial unknown vector elements. It is then convoluted with the degradation function 𝑝𝑖𝑗 – for a gamma-ray spectrum for instance this function would correspond to the resolution of the detector employed to acquire it – yielding the spectrum 𝑐𝑖 𝑝𝑖𝑗 . Each element of the vector 𝑐𝑖 is therefore spread out to its neighbor companions, losing part of its own original contents and gaining a variable contribution of them, depending upon the distance between i and j. A vector 𝑢0𝑗 – where the upper index refers to the iteration number – is assigned as an initial guess and undergoes a convolution with the same degradation function 𝑝𝑖𝑗 yielding the vector 𝑔𝑖𝑡 (left). A correcting vector 𝑓𝑗 is then obtained by a proper sum along i of the ratio 𝑐𝑖 𝑝𝑖𝑗 ∕𝑔𝑖𝑡 as sketched in Fig. 8 which is fairly self-explanatory. The arrow indicates that the LSF, with epicenter at j, scans the whole spectrum along i, computing at each j a correcting element of the vector 𝑓𝑗 . Both vectors 𝑐𝑖 𝑝𝑖𝑗 and 𝑔𝑖𝑡 have been produced by a convolution with the same degradation function. The only difference is that the last one arises from a guess vector while the first one comes from an actually observed spectrum. One can then straightforward realize that an underestimated guess for 𝑔𝑖𝑡 produces a high valued 𝑓𝑗 which corrects 𝑢𝑡𝑗 increasing it to 𝑢𝑡+1 and vice-versa. 𝑗 For images, the vectors are replaced with matrices, and the 𝑝𝑖𝑗 function LSF is replaced with PSF. Its FWHM is assigned as the spatial resolution of the acquisition system, i.e., the width of the function which degraded the primordial unknown matrix, and is used as an unfolding function to recover the primordial image. As an inverse problem, this is
Fig. 5. Synthetic radiograph of a shielding blade and its related Edge Response Function – ERF. Its derivative yields the Line Spread Function – LSF fitted with a Gaussian which FWHM is assigned as the PSF width.
as: 𝑞(𝑖) =
𝑁 1 ∑ 𝑝(𝑖, 𝑗) 𝑁 𝑗=1
(3)
where: q(i) = ERF intensity 𝑖 = 1 to M, 𝑗 = 1 to 𝑁. The resulting ERF inherits all relevant parameters from its motherimage such as w, g, and the geometric values of L and D. After differentiation a Gaussian is fitted to the result and its FWHM assigned as the LSF width [15] as shown in Fig. 5. For a given RC beam divergence w – which, as the L/D is a characteristic parameter for each acquisition system – the PSF width s corresponds to the spatial resolution, and depends upon the objectdetector gap g as shown in Fig. 6, constructed after the generated synthetic images. For reference, besides the curve arising from the actual beam divergence of 1◦ 21′ found in Ref. [12], two neighbor curves corresponding to 1◦ 10′ and 1◦ 30′ are also included in the graph. As expected, higher divergences yield steeper curves for they produced broader penumbrae and thus poorer (higher) spatial resolutions. Once a family of synthetic curves s(w, g ) is available – which can be generated as dense as necessary – the beam divergence can be determined by a pair of coordinates [𝑔𝑥 , 𝑠𝑥 ] intercepting one of them. This pair can be obtained by a set of experimental radiographs unfolded by the RL algorithm as detailed in Section 2.4. Yet, prior to proceed, it 175
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 6. Procedure to determine the beam divergence w, using a family of synthetic curves s(w, g ) – generated as dense as necessary – and a pair of coordinates [𝑔𝑥 , 𝑠𝑥 ] intercepting one of them, defining w. The central curve corresponds to the actual divergence of the facility, expressed as the angular HWHM of its Rocking Curve. It corresponds to an L/D ratio of about 20. A scheme is included for the sake of clarity.
Fig. 7. Mechanism of the RL algorithm for one dimension. Observed 𝑐𝑖 and guess 𝑢𝑡𝑗 values are convolved by the LSF degradation function𝑝𝑖𝑗 yielding 𝑐𝑖 𝑝𝑖𝑗 and 𝑔𝑖𝑡 respectively. The summation along i of their ratio produces the vector 𝑓𝑖 which corrects the last obtained vector 𝑢𝑡𝑗 .
a difficult task which can only be partially accomplished. Yet, an overall image improvement is usually achieved. The iteration proceeds until a certain condition defined by a halt criterion is reached. Different approaches have been proposed to establish this criterion – expressed in Fig. 7 as a comparison between a parameter 𝛿 and its limiting value 𝜀 – but it is generally agreed that the ultimate decision should be performed by a visual inspection. It should be stressed that after the unfolding, the primordial vector remains unknown, for the final result 𝑢𝑡𝑗 is simply its most likely representative one. At any rate, the RL algorithm is very robust and always converges to the most likely primordial vector or matrix, disregarding the initial guess. Hence, even the observed image itself, or an utterly homogeneous one can be assigned as that.
divergence may exhibit the same PSF of another one acquired under a higher divergence if its gap is larger. This assertion can be inferred from an interception of the family of curves in Fig. 6 by a horizontal line. Hence, a pair of coordinates [g, s] is required to find the curve amidst the family s(w,g) tied to the searched divergence w. This is done by qualitative and quantitative approaches as addressed in Sections 2.4.1 and 2.4.2 respectively. 2.4.1. Qualitative evaluation of the PSF The procedure to evaluate qualitatively the PSF width is sketched in Fig. 9, which is fairly self-explained. The original radiograph is deconvolved by a set of PSFs under a constant number of iterations per PSF, and the results duly recorded as image files. A visual inspection is then performed and the PSF width which yields the best image sis assigned the best PSF width 𝑠𝑥 .
2.4. Determination of the beam divergence through Richardson–Lucy’s unfolding
2.4.2. Quantitative determination of the PSF Despite the unbeatable human eye capability to recognize patterns, a visual inspection aiming at to choose the better quality between quasisimilar images is not a reliable procedure – due to the unavoidable
The width s of the PSF degrading a radiograph is ruled by the beam divergence, defined by the HWHM w of the rocking curve and the object-detector gap g. Therefore, an image acquired under a low beam 176
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 8. Computation of the correcting vector 𝑓𝑗 in the RL algorithm. The LSF centered at j scans the whole spectrum along i for each iteration. A fair estimation of the spectrum 𝑔𝑖𝑡 is achieved when the ratio 𝑐𝑖 ∕𝑔𝑖𝑡 reaches a value close to 1.
Fig. 9. Determination of the best PSF width by the Richardson–Lucy algorithm through visual inspection. After deconvolution of the degraded radiograph with the same number of iterations for each PSF width, that yielding the best image quality is assigned as the best PSF width 𝑠𝑥 .
subjectivity – when a quantitative outcome is targeted. Therefore, a procedure based on the global contrast exhibited by an image is applied through an algorithm (Algorithm G for short) [16], briefly overviewed in Section 2.4.2.1. As schematized in Fig. 10, the determination of the beam divergence is performed by the Algorithm G, fed with the original Degraded Radiograph and a set of PSFs with increasing widths. An embedded Algorithm RL unfolds the original radiograph with a given PSF per iteration. The global contrast (G-value for short) of each resulting image reaches a maximum at the best PSF width, i.e., that yielding the best image quality, being thus assigned as 𝑠𝑥.
(5)
𝛽(𝑖, 𝑗) = 1 for 𝑢(𝑖, 𝑗) > 𝑢𝑚
(6)
𝑢𝑚 = (𝑀.𝑁)−1 .
𝑀 ∑ 𝑁 ∑
𝑢(𝑖, 𝑗)
(7)
𝑖=1 𝑗=1
where: 𝑢(𝑖, 𝑗) = Pixel value at the point (i, j). M, N = No. of columns and lines of the image matrix. As stated by Eq. (7), 𝑢𝑚 is the average pixel-value of the whole image, which is used to classify all pixels within two classes: below and above it. The sum of all pixel-values occurring above and below 𝑢𝑚 constitutes the numerator and denominator of Eq. (4) respectively defining the G-value while the bump function 𝛽 (i, j) defined by Eq. (5) and Eq. (6) cuts off the pixel-values below and above 𝑢𝑚 respectively.
2.4.2.1. Brief overview of the Algorithm G. For the sake of completeness an overview of the Algorithm G is presented. This algorithm transduces the global contrast into a single number (G-value) defined as follows: [𝑀 𝑁 ] [𝑀 𝑁 ]−1 ∑∑ ∑∑ 𝐺= 𝑢(𝑖, 𝑗).𝛽(𝑖, 𝑗) . 𝑢(𝑖, 𝑗). |1 − 𝛽(𝑖, 𝑗)| (4) 𝑖=1 𝑗=1
𝛽(𝑖, 𝑗) = 0 for 𝑢(𝑖, 𝑗) ≤ 𝑢𝑚
𝑖=1 𝑗=1
177
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 10. Determination of the best PSF width using the Richardson–Lucy algorithm embedded into the Algorithm G. The original radiograph is deconvolved with several PSF widths (1 iteration per PSF), and that yielding the highest global contrast (G-value) is assigned as the best width 𝑠𝑥 .
inspection. This work targets a ratification of robustness of the RC concept through the development and application of a novel technique to determine the beam divergence. To reach this objective, the beam divergence has been determined – after the procedure shown in Section 2.4 – under two different gaps (12.1 and 19.3 mm) in order to check the internal consistency of the technique. Indeed, the beam divergence is a characteristic parameter of the reactor port and should not be affected by the object-detector gap. As depicted in Fig. 6, a pair of coordinates [𝑔, s] would specify a point in a curve s(𝑤, g ) defining thus the value of w, i.e., the HWHM of the rocking curve expressing the beam divergence. Since the abscissa g is known, one only has to determine the ordinate s, a task which is performed in this work qualitatively (by visual inspection) and quantitatively (using the Algorithm G).
This formalism is justified as follows: after a deconvolution with a given w, many pixels in the gray region shift towards darker or brighter regions due to the image improvement arising from the penumbra reduction. Hence the histogram tails raise with a consequent lowering of the remaining profile, as its area must be conserved. A higher fraction of pixels within darker and brighter zones means a higher contrast arising from a better resolution. Therefore, the sum of all pixel-values above 𝑢𝑚 will increase while that below it will decrease, resulting into a higher G-value which is computed and stored. A chosen increment is added to w and a new deconvolution is carried out. The Gvalue increases with w, reaches a maximum, and then decreases again. This behavior follows the quality of an image undergoing deconvolution; at first there is a fair improvement until a visually optimum, followed by a degradation when the iterations proceed. The w-value yielding this maximum is assigned as the best value which should be used as the PSF width for the image deconvolution. Only 1 iteration is required to achieve a smooth and reliable G-value x PSF width curve.
3.1. Best PSF width qualitatively evaluated
2.5. Acquisition of experimental radiographs The Fig. 11 depicts the radiograph system installed at the main porthole of the Argonauta reactor in Instituto de Engenharia Nuclear, Rio de Janeiro, Brazil, including some neutronic features. The objects chosen to be inspected were plastic toothed wheels. A neutron-sensitive BAS-IP ND 2040 imaging plate, using Gd2 O3 as neutron-to-ionizing radiation converter is exposed to a 3 × 105 n. cm−2 .s−1 neutron flux during 3 min and stores the image cast by the object. After a development carried out with the 50 μm spot size laser of a FUJIFILM BAS-2500 IP-Reader the images are stored in files and transferred to the freeware ImageJ which returns with a pixel intensity matrix in tiff format. 3. Results
The best PSF width is assigned as that yielding the best image, after deconvolution of the original radiograph under several values. The results – after a deconvolution with 20 iterations – are shown in Figs. 12 and 13. It is worth to stress that while the best unfolded image for object (a) is clearly defined, that for object (b) could led to some doubts amidst their companions. Indeed, besides resolution, contrast and noise, the image should keep the object shape and proportions as fairly done by FWHM = 0.47. Its right FWHM = 0.85 mm companion, despite an apparent better resolution and contrast exhibits somewhat deformed teeth as shown in Fig. 13. It is therefore a good approach to choose a test-object with features easily recognizable in order to identify eventual deviations from their original shapes. 3.2. Best PSF width quantitatively determined The procedure to determine quantitatively the best PSF width is similar to that used in the qualitative one, but much less cumbersome, subjective and time-consuming. As its companion, it also undertakes several deconvolutions using a set of PSF widths, but only 1 iteration per width, after which the global contrast of the resulting deconvolved image is determined yielding its G-value. This value is stored and another
Two small plastic toothed wheels have been chosen as test-objects, due to their good thermal neutron attenuation properties, homogeneous teeth shapes and regular gap between them, fundamental features when the quality of an image should be evaluated or ratified by visual 178
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 11. Reactor port and geometric arrangement to acquire the neutron radiographs. The imaging plate is facing directly the paraffin block but has been displaced in the drawing for clarity. Some features of the neutron beam are also displayed.
Fig. 12. Neutron radiographs of plastic toothed wheels of 25.4 (a) and 18 mm (b) diameters. Acquired at detector gaps of 12.1 and 19.3 mm respectively and unfolded with PSF widths within the range 0.15–1.5 mm, the pairs [12.1, 0.31] and [19.3, 0.47] constitutes coordinates defining the beam divergence at the family of curves s(w, g ) shown in Fig. 15.
deconvolution (also with 1 iteration) is performed with the next PSF width and the same original radiograph. After scanning a certain PSF width range – which should naturally encompass the unknown best PSF width – the profile G x PSF width is recorded as a file, which can be retrieved and plotted. All the process runs automatically and the customer has only to seek for the maximum G of the profile. Its abscissa defines the best PSF width. The Fig. 14 shows the resulting profiles for objects (a) and (b). Since the profile (b) is somewhat flatter than (a), a zoom has been added – which amplified the ripple, smoothed by a Gaussian fit – for a better definition of the peak. This slight flatness is most likely caused by the wider penumbra cast on the detector due to the larger objectdetector gap. At any rate, the determined best PSF widths fairly agree with those visually evaluated within 3.3% and 6.8%, for object (a) and (b) respectively.
3.3. Analysis of the results The coordinate pairs [12.1, 0.31] and [19.3, 0.47] obtained at the qualitative evaluation, as well as those [12.1, 0.30] and [19.3, 0.44] at the quantitative one, are plotted in the zoomed tail of Fig. 6, as shown in Fig. 15. The remaining FWHMs 0.15, 0.85 and 1.50 mm, also tested as PSF widths, remain outside the graph limits. While the first data set lands very well in the 1◦ 21′ curve, the quantitative one deviates about 3.3% and 6.8%, for the gaps 12.1 and 19.3 respectively. It is however unlikely that an intrinsically subjective visual inspection could exhibit a lower uncertainty than a quantitative technique employing an algorithm which had been earlier tested with several synthetic and experimental neutron radiographs [16]. Rather, it would be more conservative to allow an uncertainty of circa 7% for the quantitative procedure. At any rate, this uncertainty could be fairly 179
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 13. Criterion to select the best image quality. Besides resolution, contrast and noise, the image should keep the object shape and proportions as fairly done by FWHM=0.47. Its right 0.85 mm companion, despite an apparent better resolution and contrast, exhibits somewhat deformed teeth.
Fig. 14. Best PSF widths quantitatively determined for objects (a) and (b), values that fairly agree with those qualitatively evaluated within 3.3% and 6.8% respectively. Due to its flatter peak, a zoom has been added to the right graph.
acceptable taking into account the extent of error propagation occurring along the procedures required to reach the final outcome. 4. Conclusion A novel technique is presented to determine the divergence of a neutron beam aiming at the ratification of the Rocking Curve – RC concept, as recently posited by two different works. Due to their mutual agreement, soundness and consistence, it was worthwhile to cross-check the robustness of the concept through a different technique. This is performed through qualitative and quantitative procedures. In both cases, an ad hoc written Fortran 90 program generates a family of curves s(w, g ) where s is the spatial resolution, g the object-detector gap, and w the beam divergence. The experimental radiograph of an object acquired under a known gap g is then unfolded with the Richardson– Lucy algorithm using several spatial resolutions s. For the qualitative evaluation, a chosen constant number of iterations per deconvolution is carried out for each s and that 𝑠𝑥 producing the best image – visually chosen – together with g constitute the point of coordinates [g, 𝑠x ] which is plotted in the graph containing the family of curves s(w, g ) hitting one of them. Since each curve is tied to a given beam divergence, the
Fig. 15. Zoomed tail of Fig. 6. Radiographs of objects (a) and (b), as seen in Fig. 12, taken at gaps 12.1 and 19.3 mm respectively have been unfolded with several PSF widths. The best images – after a visual evaluation – occur at 0.31 and 0.47 mm, hitting the 1◦ 21′ curve, while quantitative results occur at 0.30 and 0.44 mm.
180
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
hit one defines the divergence w at which the original radiograph had been acquired. The quantitative procedure employs a similar technique except that 𝑠x is determined through the maximum global contrast exhibited by the images deconvolved with only 1 iteration. The beam divergence determined in this work agrees within 7% with those formerly measured – using utterly different techniques – an outcome which indicates the soundness of the presented method and points towards the robustness of the rocking curve approach.
[7] L. Greim, M. Greim, H.W. Schmitz, G.W. Schumacher, Computer controlled microdensitometer and some applications, neutron radiography, in: Proceedings of the 2nd World Conference, June 16–20, Paris, France, D. Reydel Publishing Company, Dordrecht, 1986, pp. 669–677. [8] H. Kobayashi, Neutron radiography research in Rikkyo reactor, in: Proceedings of the 1st Asian Symposium on Research Reactors, November 18–21, Rikkyo University, Tokio, Japan, 1986, pp. 360–370. [9] H. Kobayashi, H. Wakao, Accurate measurement of l, d and l/d for divergent collimators, in: Proceedings of the 3rd World Conference, May 14–18, 1989 Osaka, Japan, Kulwer Academic Publishers, Dordrecht, 1989, pp. 885–892. [10] Y. Nir-El, W.L. Whittemore, Accurate measurement of L, D and L/D for divergent collimators, in: Proceedings of the 3rd World Conference, May 14–18, 1989 Osaka, Japan, Kulwer Academic Publishers, Dordrecht, 1989, pp. 885–892. [11] G.L. Almeida, M.I. Silvani, R.C. Furieri, M.J. Goncalves, R.T. Lopes, Braz. J. Phys. 35 (3B) (2005) 771–774. [12] E.S. Souza, G.L. Almeida, R.T. Lopes, Nucl. Instrum. Methods Phys. Res. A 838 (2016) 129–136. [13] M.I.S. Souza, Computer-Aided Tomography with Thermal Neutrons and a Proportional Gaseous Position-Sensitive Detector (D.Sc. thesis), Universidade Federal do Rio de Janeiro - UFRJ, COPPE, 2001 (in Portuguese). [14] S.W. Smith, The Scientist and Engineer’s Guide to Digital Signal Processing, California Technical Publishing, 2011. [15] ASTM E 1441-95 and 1570-95a, Non-destructive Testing, Radiation Methods, Computed Tomography, Guide for Imaging and Practice Examination, ISO/TC 135/SC5, N118, 1996. [16] G.L. Almeida, M.I. Souza, A novel algorithm for blind deconvolution applied to the improvement of radiographic images, AIP Conf. Proc. 1529 (2013) 95. http: //dx.doi.org/10.1063/1.4804094.
Acknowledgments This work was partially supported by Brazilian agencies Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). References [1] W.H. Richardson, J. Opt. Soc. Amer. 62 (1) (1972) 55–59. [2] B.L. Lucy, Astron. J. 79 (6) (1974) 745–754. [3] H. Berger, Neutron Radiography. Methods, Capabilities, and Applications, Elsevier Elsevier Publishing Co., Amsterdam, 1965. [4] J.D. Domanus, Practical Neutron Nadiography, Kulwer Academic Publishers, Dordrecht, 1992. [5] J.P. Barton, Materials Evaluation (1967) 45A. [6] ASTM E803-81, Standard method for determining the L/D ratio of neutron radiographic beams.
181
Contents lists available at ScienceDirect
Nuclear Inst. and Methods in Physics Research, A journal homepage: www.elsevier.com/locate/nima
Determination of a neutron beam divergence after the rocking curve concept using Richardson–Lucy’s unfolding algorithm E.S. Souza a ,∗, M.I. Silvani b , G.L. Almeida b , R.T. Lopes a a
Universidade Federal do Rio de Janeiro, COPPE, Centro de Tecnologia, Cidade Universitaria Bloco G, Ilha do Fundao, 21945-970 Rio de Janeiro - RJ, Brazil Instituto de Engenharia Nuclear, Reator Argonauta - CNEN, Rua Helio de Almeida 75, Cidade Universitária, Ilha do Fundao Caixa Postal 68550 CEP 21941-972 Rio de Janeiro - RJ, Brazil b
ARTICLE
INFO
Keywords: Neutron beam divergence Rocking curve L/D Richardson–Lucy
ABSTRACT Thermal neutron radiographs acquired under high beam divergences suffer the impact of intense penumbrae degrading their final quality. As a divergence reduction is not always feasible, one possible alternative is an a posteriori treatment to restore the degraded images, such as the Richardson–Lucy — RL unfolding algorithm. Such a procedure requires the characterization of the spoiling agent Point Spread function – PSF in order to apply it in the inverse way. It can be deduced from the blur in the radiograph of a shielding blade edge. This blur depends upon the beam divergence and the object–detector gap. Due to the complex scattering of neutrons along a reactor channel, it is usual to express this divergence as the inverse of an L/D ratio. A novel approach based on the concept of Rocking Curve – RC, a term borrowed from the X-ray diffraction field, has been recently proposed which yielded slightly better quantitative results. After this concept, every sub-element of a surface source emits neutrons anisotropically following a bell-shaped profile. The RC angular semi-width incorporating neutron scattering and geometrical blur, is assigned as the beam divergence. The present work aims at its assessment through a quantitative determination ratified by a qualitative evaluation of radiographs unfolded by the RL algorithm. Its main purpose is an additional ratification of the soundness, consistency and robustness of the RC concept by comparison with formerly obtained quantitative results. In spite of the utterly different approaches and techniques, the outcome has corroborated the novel concept. All data treatment is simple and performed by an ad hoc written Fortran 90 program embedding the required algorithms.
1. Introduction
depends upon the beam divergence w, being amplified by the objectdetector gap g.
Radiographic images acquired with thermal neutrons are infested by plagues such as beam divergence, poor detector resolution, neutron scattering, statistical fluctuation and electronic noise which degrade their final quality. All of them are impossible to eliminate or at least hard to mitigate at their origin due to their fundamental features or technological and engineering constraints. Yet, although not possible to restore the primordial image, it is feasible to improve them through a post-acquisition treatment such as an unfolding with the Richardson– Lucy — RL algorithm [1,2]. This treatment requires however the characterization of the spoiling agent Point Spread Function–PSF in order to apply it in the inverse way. The PSF is the response of an image acquisition system to a pointlike input. For isotropic systems it has an intensity distribution with a radial symmetry engulfing all the spoiling agents. Although hard to characterize, it can be fairly expressed by the rotation of a Gaussian around its epicenter, and its standard deviation deduced from the blur cast on a detector by the straight edge of shielding blade. This blur
As for the beam divergence itself, it is worthwhile to define its meaning within the scope of this work as follows. A flat surface source emitting only in the perpendicular direction, should not produce any penumbra at the detector, independently of the object-detector gap or its thickness. As early stated by Berger [3] – quoted in Domanus [4] – ‘‘. . . the present state of art is such that parallel beams are definitely preferred’’. This assertion was later on refuted by Barton [5] who concluded that a divergent collimator produces better images. This disagreement seems to be a matter of misinterpretation. Perhaps Berger had in mind an ideal flat source as above described and not a real one where each source sub-element would emit neutrons isotropically. Or perhaps he was actually referring to quasi-parallel beams, as those tied to high L/D ratios. Under such a circumstance, a divergent collimator would naturally produce a better image, as stated by Barton, because in this case, the ratio of the source–detector clearance L to source size D
∗ Corresponding author. E-mail addresses: [email protected] (E.S. Souza), [email protected] (M.I. Silvani), [email protected] (G.L. Almeida), [email protected] (R.T. Lopes).
https://doi.org/10.1016/j.nima.2018.08.055 Received 15 May 2018; Received in revised form 2 August 2018; Accepted 18 August 2018 Available online xxxx 0168-9002/© 2018 Published by Elsevier B.V.
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
would increase towards that exhibited by an ideal point source, where the L/D ratio would obviously raise to infinite. An ideal point-source as well, would not blur the image, and thus is treated in this work as nondivergent, in spite of its conical beam. In a real neutron acquisition system however – where a virtual flat source could be visualized somewhere in the reactor neutron channel – scattering causes some neutrons to hit the detector at different angles, as if they were delivered from several sources at different distances and intensities. The overall impact on the quality of final image is equivalent to that caused by a single source of unknown size D located at an as well unknown distance L from the detector. As comprehensively addressed by Domanus [4], other techniques e.g., [6–10] rather than simple length measurements should be carried out to measure the effective L/D, which is always lower than the purely geometric one and rules the quality of the final image. Hence, disregarding the specific geometric arrangement at any facility, the images would exhibit equivalent qualities if they were acquired under the same effective L/D. The approach condensing the spoiling agents into a single L/D parameter overcomes the complex task of dealing with them, mainly with neutron scattering. Recent works have been proposed [11,12] positing the concept of Rocking Curve – RC, a designation borrowed from the X-ray diffraction field. It assumes that the source emit neutrons anisotropically with intensities ruled by a 3D-Gaussian distribution, and the beam divergence expressed by its angular half width at half maximum – HWHM designated as w in this work. For the sake of completeness, some features of those works including drawings, sketches and equations are shared by the present work. Although both earlier works share the same concept, they employ utterly different approaches and techniques. While Ref. [11] employs a position sensitive detector, a slit collimator and data treatment involving unfolding and extrapolation, Ref. [12] requires the neutron radiograph of a shielding blade edge, and an ad hoc written algorithm. Despite these differences, the results agree within 7% as follows: Reference RC HWHM (min)
[11] 75.91 ± 1.52
the searched one tied to the RC HWHM 𝑤0 . Further details are presented in Sections 2.4 and 2.5. An illustrative concept of the RC is sketched in Fig. 1, where a high divergence corresponds to a broad Gaussian (large w) while a zero divergence would be expressed by an infinitely narrow one, emitting neutrons solely in the axial direction. 2. Materials and methods The direct experimental determination of the PSF is very hard due to the low counting statistics imposed by its small size. Instead, it is replaced by the measurement of the Line Spread function – LSF, the response of a system to a line which can be obtained by differentiation of the Edge Response Function – ERF. This function is provided by the radiograph of a shielding blade edge. A Gaussian is then fitted to the LSF distribution and its FWHM assigned as the PSF width, since LSF and PSF share the same FWHM when the distribution follows a Gaussian profile [14] as depicted in Fig. 2. To overcome the related radiological burden, reactor operation costs and a cumbersome work, the family of curves s(w, g ) is obtained by an ad hoc written algorithm which generates synthetic radiographs and carries out the required treatments. 2.1. Generation of the synthetic images To achieve this task a virtual 2D detector emulating a real one such as an imaging plate is positioned at the end of a neutron channel perpendicularly to its axis, as sketched in Fig. 3. Somewhere in the neutron channel, a flat neutron source filling its whole cross-section would emit neutrons which could or not be intercepted by a shielding blade placed at a chosen distance g from the detector. Since its edge is aligned with the vertical axis dividing the image into two regions, the left region would be shielded while its right companion would be hit by the neutrons coming from the source. This ideal situation would occur only if the beam divergence or g were zero. For all other cases the penumbra would invade the neighborhood of the border line blurring the otherwise sharp edge. The intensity of the neutrons hitting the detector is assumed to depend upon the angle 𝛷 of their paths to the normal to the source as illustrated by the pictured bell-shaped surface. Once a spatial detector resolution 𝛿 is defined, an MxN matrix is assigned to the image, so that 𝛿.M and 𝛿.N represent its width and height respectively. Each pixel intensity is obtained by the summation of all neutrons coming from the source subdivided like the detector into KxH square elements of size 𝛼. So, 𝛼.K and 𝛼.H represent the width and height of the neutron channel and its virtual source. For a real source-detector set, the pixel intensity would depend upon the source intensity, detection efficiency and exposure time, but since that – for tiff images – it should not surpass the limit of 65,535, after summation, the pixel intensities as defined by Eq. (1) are normalized to this limit.
[12] 81.20 ± 1.18
Besides this fair agreement, the single blade required by Ref. [12], instead of an elaborated test-object to measure the L/D, makes the method affordable and attractive due to its simplicity. Furthermore, its basic algorithm used to measure the RC has been also used to measure the L/D agreeing fairly with Ref. [13] – based on neutron flux measurements – and exhibiting a lower uncertainty than the standard method presented in Ref. [6] as follows: Reference L/D ratio L/D uncertainty %
[6] 12.5–16.7
[12] 20.14 ± 1.17 5.8
[13] 19.2 ± 0.71 3.7
The above data refer to the neutron port of the Argonauta research reactor at Instituto de Engenharia Nuclear – CNEN, Rio de Janeiro Brazil, where the measurements [11–13] have been done. In spite of the better uncertainty exhibited by Ref. [13], it is less employed due to its timeconsuming. Taking into account this performance, and considering that the RC concept has not been yet comprehensively and exhaustively crosschecked, it is an advisable policy to ratify its soundness, consistency and robustness by another approach. Within this frame, the present work addresses the determination of the beam divergence through unfolding of a neutron radiograph with the RL algorithm. This is done by an analysis of radiographs unfolded with different PSF widths s, which selects the better amidst them. Its related PSF width 𝑠𝑥 and the object-detector gap 𝑔0 used to get the original radiograph constitute the coordinates [g0 , 𝑠x ] from which the related beam divergence 𝑤0 will be inferred. To accomplish this task, a family of synthetic curves s(w, g ) is generated and the curve s(𝑤0 , g) hit by [g0 , 𝑠x ] is assigned as
𝑝(𝑖, 𝑗) =
𝐾 ∑
𝛺(𝑖, 𝑘).
𝑘=1
𝛷𝑤 = 𝑤.][𝑙𝑛4
𝐻 ∑
{ [ ]2 } exp −0.5 𝛷 (𝑖, 𝑗, 𝑘, ℎ) ∕𝛷𝑤
(1)
ℎ=1
(2)
where: p(𝑖, 𝑗) = Pixel intensity on the detector: 𝑖 = 1 to M, 𝑗 = 1 to 𝑁. 𝛺(𝑖, 𝑘) = Bump function: 0 if the neutron hits the shielding, 1 otherwise. 𝛷(𝑖, 𝑗, 𝑘, ℎ) = Angle between the normal to the source and the straight line connecting the points (𝑖, 𝑗) on the detector and (𝑘, ℎ) on the source. 𝛷𝑤 = Angular standard deviation of the generatrix Gaussian, for the 3D Rocking Curve. w = Angular Half-Width at Half Maximum − HWHM of the Gaussian. K = No. of elements along the source width. H = No. of elements along the source height. M = No. of pixels along the image width on the detector. N = No. of pixels along the image height on the detector. 173
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 1. Concept of Rocking Curve as an anisotropic emission characterizing the beam divergence defined by its half width at half maximum – HWHM w.
Fig. 2. Response of an image acquisition system to a point and a line. Since the Point Spread Function – PSF is hard to measure experimentally it is replaced by the Line Spread Function – LSF. For a Gaussian distribution, they share the same FWHM. Fig. 3. Generation of synthetic images. Each pixel on the detector – one is highlighted – collects and integrates the radiation coming from all source elements. Unless blocked by the shielding blade, its intensity depends upon the angle 𝛷 of their paths to the normal to the source as illustrated by the pictured bell-shaped surface.
It should be pointed out that 𝛺(𝑖, 𝑘) does not depend upon j and h, because the shielding edge has a vertical orientation. As sketched in Fig. 4, which is fairly self-explanatory, the bell-shaped RC – represented by the rotation of a Gaussian around its axis – rules the flux of neutron emitted by each source element. Unlike an isotropic emission, the flux intensity depends upon the angle 𝛷(𝑖, 𝑗, 𝑘, ℎ) between the normal to the source and the straight line connecting the pixel (𝑖, 𝑗) on the detector and the element (𝑘, ℎ) on the source. Source and detector share the same dimensions of the neutron channel, but different pixel sizes. Both can be freely chosen as input data, but in order to emulate the detector actually used to acquire the experimental radiographs, a 0.05 mm resolution has been assigned to the virtual one. As for the source, the elements can be set as small as wished, provided that the choice does not enlarge the processing time to unreasonable values. Other parameters which should be entered as input data for the 3D analytical geometry treatment are w, g, L, and D.
As it is assumed that the RC can be represented by a Gaussian function, some limit should be set to discard insignificant contributions arising from its tails. However to keep a conservative approach, all contributions coming from angles below F.w are taken into account. An F -value of 5.0 has been employed in this work. 2.2. Determination of the ERF, LSF and PSF width for the synthetic images An ERF is obtained by averaging the pixel intensities along a line parallel to the edge shadow, both for synthetic and experimental images. This averaging encompassing several pixels along a row aims at to smooth the unavoidable noise and statistical fluctuation, and is written 174
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 4. 3D Rocking Curve as an extended concept borrowed from the X-ray diffraction field using a Gaussian as generatrix. The flux intensities 𝑆𝑎 and 𝑆𝑏 are tied to the respective angles 𝛷𝑎 and 𝛷𝑏 .
is worth to review briefly the mechanism of this algorithm as depicted in Figs. 7 and 8 due to its paramount role in this work. 2.3. Brief overview of the Richardson–Lucy unfolding algorithm In order to simplify the description it is convenient to choose a 1D example, such as a gamma-ray spectrum. In the absence of noise of any kind, and with an ideal zero-resolution detector, the gamma-rays would appear as vertical lines along the energy axis. If these conditions are not fulfilled, the contents (counts/channel) of every channel leaks to neighbor companions contaminating them with contribution of alien energies. In consequence of this convolution, the otherwise vertical lines are degraded to the usual peak shape. Yet, the information is not lost – just concealed – being capable to be retrieved by a proper treatment as schematized in Fig. 7. The RL algorithm employs an iterative procedure to determine the 𝑡 vector 𝑢𝑡+1 𝑗 from its predecessor 𝑢𝑗 through the correcting vector 𝑓𝑗 . Hence, when 𝑓𝑗 approaches the unitary vector, the spectra become more similar to each other after an asymptotical behavior. Expressing in words, what is already precisely described by the equations, the observed spectrum 𝑐𝑖 (right) is considered as if it were a non-degraded one. The vertical strips represent the primordial unknown vector elements. It is then convoluted with the degradation function 𝑝𝑖𝑗 – for a gamma-ray spectrum for instance this function would correspond to the resolution of the detector employed to acquire it – yielding the spectrum 𝑐𝑖 𝑝𝑖𝑗 . Each element of the vector 𝑐𝑖 is therefore spread out to its neighbor companions, losing part of its own original contents and gaining a variable contribution of them, depending upon the distance between i and j. A vector 𝑢0𝑗 – where the upper index refers to the iteration number – is assigned as an initial guess and undergoes a convolution with the same degradation function 𝑝𝑖𝑗 yielding the vector 𝑔𝑖𝑡 (left). A correcting vector 𝑓𝑗 is then obtained by a proper sum along i of the ratio 𝑐𝑖 𝑝𝑖𝑗 ∕𝑔𝑖𝑡 as sketched in Fig. 8 which is fairly self-explanatory. The arrow indicates that the LSF, with epicenter at j, scans the whole spectrum along i, computing at each j a correcting element of the vector 𝑓𝑗 . Both vectors 𝑐𝑖 𝑝𝑖𝑗 and 𝑔𝑖𝑡 have been produced by a convolution with the same degradation function. The only difference is that the last one arises from a guess vector while the first one comes from an actually observed spectrum. One can then straightforward realize that an underestimated guess for 𝑔𝑖𝑡 produces a high valued 𝑓𝑗 which corrects 𝑢𝑡𝑗 increasing it to 𝑢𝑡+1 and vice-versa. 𝑗 For images, the vectors are replaced with matrices, and the 𝑝𝑖𝑗 function LSF is replaced with PSF. Its FWHM is assigned as the spatial resolution of the acquisition system, i.e., the width of the function which degraded the primordial unknown matrix, and is used as an unfolding function to recover the primordial image. As an inverse problem, this is
Fig. 5. Synthetic radiograph of a shielding blade and its related Edge Response Function – ERF. Its derivative yields the Line Spread Function – LSF fitted with a Gaussian which FWHM is assigned as the PSF width.
as: 𝑞(𝑖) =
𝑁 1 ∑ 𝑝(𝑖, 𝑗) 𝑁 𝑗=1
(3)
where: q(i) = ERF intensity 𝑖 = 1 to M, 𝑗 = 1 to 𝑁. The resulting ERF inherits all relevant parameters from its motherimage such as w, g, and the geometric values of L and D. After differentiation a Gaussian is fitted to the result and its FWHM assigned as the LSF width [15] as shown in Fig. 5. For a given RC beam divergence w – which, as the L/D is a characteristic parameter for each acquisition system – the PSF width s corresponds to the spatial resolution, and depends upon the objectdetector gap g as shown in Fig. 6, constructed after the generated synthetic images. For reference, besides the curve arising from the actual beam divergence of 1◦ 21′ found in Ref. [12], two neighbor curves corresponding to 1◦ 10′ and 1◦ 30′ are also included in the graph. As expected, higher divergences yield steeper curves for they produced broader penumbrae and thus poorer (higher) spatial resolutions. Once a family of synthetic curves s(w, g ) is available – which can be generated as dense as necessary – the beam divergence can be determined by a pair of coordinates [𝑔𝑥 , 𝑠𝑥 ] intercepting one of them. This pair can be obtained by a set of experimental radiographs unfolded by the RL algorithm as detailed in Section 2.4. Yet, prior to proceed, it 175
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 6. Procedure to determine the beam divergence w, using a family of synthetic curves s(w, g ) – generated as dense as necessary – and a pair of coordinates [𝑔𝑥 , 𝑠𝑥 ] intercepting one of them, defining w. The central curve corresponds to the actual divergence of the facility, expressed as the angular HWHM of its Rocking Curve. It corresponds to an L/D ratio of about 20. A scheme is included for the sake of clarity.
Fig. 7. Mechanism of the RL algorithm for one dimension. Observed 𝑐𝑖 and guess 𝑢𝑡𝑗 values are convolved by the LSF degradation function𝑝𝑖𝑗 yielding 𝑐𝑖 𝑝𝑖𝑗 and 𝑔𝑖𝑡 respectively. The summation along i of their ratio produces the vector 𝑓𝑖 which corrects the last obtained vector 𝑢𝑡𝑗 .
a difficult task which can only be partially accomplished. Yet, an overall image improvement is usually achieved. The iteration proceeds until a certain condition defined by a halt criterion is reached. Different approaches have been proposed to establish this criterion – expressed in Fig. 7 as a comparison between a parameter 𝛿 and its limiting value 𝜀 – but it is generally agreed that the ultimate decision should be performed by a visual inspection. It should be stressed that after the unfolding, the primordial vector remains unknown, for the final result 𝑢𝑡𝑗 is simply its most likely representative one. At any rate, the RL algorithm is very robust and always converges to the most likely primordial vector or matrix, disregarding the initial guess. Hence, even the observed image itself, or an utterly homogeneous one can be assigned as that.
divergence may exhibit the same PSF of another one acquired under a higher divergence if its gap is larger. This assertion can be inferred from an interception of the family of curves in Fig. 6 by a horizontal line. Hence, a pair of coordinates [g, s] is required to find the curve amidst the family s(w,g) tied to the searched divergence w. This is done by qualitative and quantitative approaches as addressed in Sections 2.4.1 and 2.4.2 respectively. 2.4.1. Qualitative evaluation of the PSF The procedure to evaluate qualitatively the PSF width is sketched in Fig. 9, which is fairly self-explained. The original radiograph is deconvolved by a set of PSFs under a constant number of iterations per PSF, and the results duly recorded as image files. A visual inspection is then performed and the PSF width which yields the best image sis assigned the best PSF width 𝑠𝑥 .
2.4. Determination of the beam divergence through Richardson–Lucy’s unfolding
2.4.2. Quantitative determination of the PSF Despite the unbeatable human eye capability to recognize patterns, a visual inspection aiming at to choose the better quality between quasisimilar images is not a reliable procedure – due to the unavoidable
The width s of the PSF degrading a radiograph is ruled by the beam divergence, defined by the HWHM w of the rocking curve and the object-detector gap g. Therefore, an image acquired under a low beam 176
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 8. Computation of the correcting vector 𝑓𝑗 in the RL algorithm. The LSF centered at j scans the whole spectrum along i for each iteration. A fair estimation of the spectrum 𝑔𝑖𝑡 is achieved when the ratio 𝑐𝑖 ∕𝑔𝑖𝑡 reaches a value close to 1.
Fig. 9. Determination of the best PSF width by the Richardson–Lucy algorithm through visual inspection. After deconvolution of the degraded radiograph with the same number of iterations for each PSF width, that yielding the best image quality is assigned as the best PSF width 𝑠𝑥 .
subjectivity – when a quantitative outcome is targeted. Therefore, a procedure based on the global contrast exhibited by an image is applied through an algorithm (Algorithm G for short) [16], briefly overviewed in Section 2.4.2.1. As schematized in Fig. 10, the determination of the beam divergence is performed by the Algorithm G, fed with the original Degraded Radiograph and a set of PSFs with increasing widths. An embedded Algorithm RL unfolds the original radiograph with a given PSF per iteration. The global contrast (G-value for short) of each resulting image reaches a maximum at the best PSF width, i.e., that yielding the best image quality, being thus assigned as 𝑠𝑥.
(5)
𝛽(𝑖, 𝑗) = 1 for 𝑢(𝑖, 𝑗) > 𝑢𝑚
(6)
𝑢𝑚 = (𝑀.𝑁)−1 .
𝑀 ∑ 𝑁 ∑
𝑢(𝑖, 𝑗)
(7)
𝑖=1 𝑗=1
where: 𝑢(𝑖, 𝑗) = Pixel value at the point (i, j). M, N = No. of columns and lines of the image matrix. As stated by Eq. (7), 𝑢𝑚 is the average pixel-value of the whole image, which is used to classify all pixels within two classes: below and above it. The sum of all pixel-values occurring above and below 𝑢𝑚 constitutes the numerator and denominator of Eq. (4) respectively defining the G-value while the bump function 𝛽 (i, j) defined by Eq. (5) and Eq. (6) cuts off the pixel-values below and above 𝑢𝑚 respectively.
2.4.2.1. Brief overview of the Algorithm G. For the sake of completeness an overview of the Algorithm G is presented. This algorithm transduces the global contrast into a single number (G-value) defined as follows: [𝑀 𝑁 ] [𝑀 𝑁 ]−1 ∑∑ ∑∑ 𝐺= 𝑢(𝑖, 𝑗).𝛽(𝑖, 𝑗) . 𝑢(𝑖, 𝑗). |1 − 𝛽(𝑖, 𝑗)| (4) 𝑖=1 𝑗=1
𝛽(𝑖, 𝑗) = 0 for 𝑢(𝑖, 𝑗) ≤ 𝑢𝑚
𝑖=1 𝑗=1
177
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 10. Determination of the best PSF width using the Richardson–Lucy algorithm embedded into the Algorithm G. The original radiograph is deconvolved with several PSF widths (1 iteration per PSF), and that yielding the highest global contrast (G-value) is assigned as the best width 𝑠𝑥 .
inspection. This work targets a ratification of robustness of the RC concept through the development and application of a novel technique to determine the beam divergence. To reach this objective, the beam divergence has been determined – after the procedure shown in Section 2.4 – under two different gaps (12.1 and 19.3 mm) in order to check the internal consistency of the technique. Indeed, the beam divergence is a characteristic parameter of the reactor port and should not be affected by the object-detector gap. As depicted in Fig. 6, a pair of coordinates [𝑔, s] would specify a point in a curve s(𝑤, g ) defining thus the value of w, i.e., the HWHM of the rocking curve expressing the beam divergence. Since the abscissa g is known, one only has to determine the ordinate s, a task which is performed in this work qualitatively (by visual inspection) and quantitatively (using the Algorithm G).
This formalism is justified as follows: after a deconvolution with a given w, many pixels in the gray region shift towards darker or brighter regions due to the image improvement arising from the penumbra reduction. Hence the histogram tails raise with a consequent lowering of the remaining profile, as its area must be conserved. A higher fraction of pixels within darker and brighter zones means a higher contrast arising from a better resolution. Therefore, the sum of all pixel-values above 𝑢𝑚 will increase while that below it will decrease, resulting into a higher G-value which is computed and stored. A chosen increment is added to w and a new deconvolution is carried out. The Gvalue increases with w, reaches a maximum, and then decreases again. This behavior follows the quality of an image undergoing deconvolution; at first there is a fair improvement until a visually optimum, followed by a degradation when the iterations proceed. The w-value yielding this maximum is assigned as the best value which should be used as the PSF width for the image deconvolution. Only 1 iteration is required to achieve a smooth and reliable G-value x PSF width curve.
3.1. Best PSF width qualitatively evaluated
2.5. Acquisition of experimental radiographs The Fig. 11 depicts the radiograph system installed at the main porthole of the Argonauta reactor in Instituto de Engenharia Nuclear, Rio de Janeiro, Brazil, including some neutronic features. The objects chosen to be inspected were plastic toothed wheels. A neutron-sensitive BAS-IP ND 2040 imaging plate, using Gd2 O3 as neutron-to-ionizing radiation converter is exposed to a 3 × 105 n. cm−2 .s−1 neutron flux during 3 min and stores the image cast by the object. After a development carried out with the 50 μm spot size laser of a FUJIFILM BAS-2500 IP-Reader the images are stored in files and transferred to the freeware ImageJ which returns with a pixel intensity matrix in tiff format. 3. Results
The best PSF width is assigned as that yielding the best image, after deconvolution of the original radiograph under several values. The results – after a deconvolution with 20 iterations – are shown in Figs. 12 and 13. It is worth to stress that while the best unfolded image for object (a) is clearly defined, that for object (b) could led to some doubts amidst their companions. Indeed, besides resolution, contrast and noise, the image should keep the object shape and proportions as fairly done by FWHM = 0.47. Its right FWHM = 0.85 mm companion, despite an apparent better resolution and contrast exhibits somewhat deformed teeth as shown in Fig. 13. It is therefore a good approach to choose a test-object with features easily recognizable in order to identify eventual deviations from their original shapes. 3.2. Best PSF width quantitatively determined The procedure to determine quantitatively the best PSF width is similar to that used in the qualitative one, but much less cumbersome, subjective and time-consuming. As its companion, it also undertakes several deconvolutions using a set of PSF widths, but only 1 iteration per width, after which the global contrast of the resulting deconvolved image is determined yielding its G-value. This value is stored and another
Two small plastic toothed wheels have been chosen as test-objects, due to their good thermal neutron attenuation properties, homogeneous teeth shapes and regular gap between them, fundamental features when the quality of an image should be evaluated or ratified by visual 178
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 11. Reactor port and geometric arrangement to acquire the neutron radiographs. The imaging plate is facing directly the paraffin block but has been displaced in the drawing for clarity. Some features of the neutron beam are also displayed.
Fig. 12. Neutron radiographs of plastic toothed wheels of 25.4 (a) and 18 mm (b) diameters. Acquired at detector gaps of 12.1 and 19.3 mm respectively and unfolded with PSF widths within the range 0.15–1.5 mm, the pairs [12.1, 0.31] and [19.3, 0.47] constitutes coordinates defining the beam divergence at the family of curves s(w, g ) shown in Fig. 15.
deconvolution (also with 1 iteration) is performed with the next PSF width and the same original radiograph. After scanning a certain PSF width range – which should naturally encompass the unknown best PSF width – the profile G x PSF width is recorded as a file, which can be retrieved and plotted. All the process runs automatically and the customer has only to seek for the maximum G of the profile. Its abscissa defines the best PSF width. The Fig. 14 shows the resulting profiles for objects (a) and (b). Since the profile (b) is somewhat flatter than (a), a zoom has been added – which amplified the ripple, smoothed by a Gaussian fit – for a better definition of the peak. This slight flatness is most likely caused by the wider penumbra cast on the detector due to the larger objectdetector gap. At any rate, the determined best PSF widths fairly agree with those visually evaluated within 3.3% and 6.8%, for object (a) and (b) respectively.
3.3. Analysis of the results The coordinate pairs [12.1, 0.31] and [19.3, 0.47] obtained at the qualitative evaluation, as well as those [12.1, 0.30] and [19.3, 0.44] at the quantitative one, are plotted in the zoomed tail of Fig. 6, as shown in Fig. 15. The remaining FWHMs 0.15, 0.85 and 1.50 mm, also tested as PSF widths, remain outside the graph limits. While the first data set lands very well in the 1◦ 21′ curve, the quantitative one deviates about 3.3% and 6.8%, for the gaps 12.1 and 19.3 respectively. It is however unlikely that an intrinsically subjective visual inspection could exhibit a lower uncertainty than a quantitative technique employing an algorithm which had been earlier tested with several synthetic and experimental neutron radiographs [16]. Rather, it would be more conservative to allow an uncertainty of circa 7% for the quantitative procedure. At any rate, this uncertainty could be fairly 179
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
Fig. 13. Criterion to select the best image quality. Besides resolution, contrast and noise, the image should keep the object shape and proportions as fairly done by FWHM=0.47. Its right 0.85 mm companion, despite an apparent better resolution and contrast, exhibits somewhat deformed teeth.
Fig. 14. Best PSF widths quantitatively determined for objects (a) and (b), values that fairly agree with those qualitatively evaluated within 3.3% and 6.8% respectively. Due to its flatter peak, a zoom has been added to the right graph.
acceptable taking into account the extent of error propagation occurring along the procedures required to reach the final outcome. 4. Conclusion A novel technique is presented to determine the divergence of a neutron beam aiming at the ratification of the Rocking Curve – RC concept, as recently posited by two different works. Due to their mutual agreement, soundness and consistence, it was worthwhile to cross-check the robustness of the concept through a different technique. This is performed through qualitative and quantitative procedures. In both cases, an ad hoc written Fortran 90 program generates a family of curves s(w, g ) where s is the spatial resolution, g the object-detector gap, and w the beam divergence. The experimental radiograph of an object acquired under a known gap g is then unfolded with the Richardson– Lucy algorithm using several spatial resolutions s. For the qualitative evaluation, a chosen constant number of iterations per deconvolution is carried out for each s and that 𝑠𝑥 producing the best image – visually chosen – together with g constitute the point of coordinates [g, 𝑠x ] which is plotted in the graph containing the family of curves s(w, g ) hitting one of them. Since each curve is tied to a given beam divergence, the
Fig. 15. Zoomed tail of Fig. 6. Radiographs of objects (a) and (b), as seen in Fig. 12, taken at gaps 12.1 and 19.3 mm respectively have been unfolded with several PSF widths. The best images – after a visual evaluation – occur at 0.31 and 0.47 mm, hitting the 1◦ 21′ curve, while quantitative results occur at 0.30 and 0.44 mm.
180
E.S. Souza et al.
Nuclear Inst. and Methods in Physics Research, A 908 (2018) 172–181
hit one defines the divergence w at which the original radiograph had been acquired. The quantitative procedure employs a similar technique except that 𝑠x is determined through the maximum global contrast exhibited by the images deconvolved with only 1 iteration. The beam divergence determined in this work agrees within 7% with those formerly measured – using utterly different techniques – an outcome which indicates the soundness of the presented method and points towards the robustness of the rocking curve approach.
[7] L. Greim, M. Greim, H.W. Schmitz, G.W. Schumacher, Computer controlled microdensitometer and some applications, neutron radiography, in: Proceedings of the 2nd World Conference, June 16–20, Paris, France, D. Reydel Publishing Company, Dordrecht, 1986, pp. 669–677. [8] H. Kobayashi, Neutron radiography research in Rikkyo reactor, in: Proceedings of the 1st Asian Symposium on Research Reactors, November 18–21, Rikkyo University, Tokio, Japan, 1986, pp. 360–370. [9] H. Kobayashi, H. Wakao, Accurate measurement of l, d and l/d for divergent collimators, in: Proceedings of the 3rd World Conference, May 14–18, 1989 Osaka, Japan, Kulwer Academic Publishers, Dordrecht, 1989, pp. 885–892. [10] Y. Nir-El, W.L. Whittemore, Accurate measurement of L, D and L/D for divergent collimators, in: Proceedings of the 3rd World Conference, May 14–18, 1989 Osaka, Japan, Kulwer Academic Publishers, Dordrecht, 1989, pp. 885–892. [11] G.L. Almeida, M.I. Silvani, R.C. Furieri, M.J. Goncalves, R.T. Lopes, Braz. J. Phys. 35 (3B) (2005) 771–774. [12] E.S. Souza, G.L. Almeida, R.T. Lopes, Nucl. Instrum. Methods Phys. Res. A 838 (2016) 129–136. [13] M.I.S. Souza, Computer-Aided Tomography with Thermal Neutrons and a Proportional Gaseous Position-Sensitive Detector (D.Sc. thesis), Universidade Federal do Rio de Janeiro - UFRJ, COPPE, 2001 (in Portuguese). [14] S.W. Smith, The Scientist and Engineer’s Guide to Digital Signal Processing, California Technical Publishing, 2011. [15] ASTM E 1441-95 and 1570-95a, Non-destructive Testing, Radiation Methods, Computed Tomography, Guide for Imaging and Practice Examination, ISO/TC 135/SC5, N118, 1996. [16] G.L. Almeida, M.I. Souza, A novel algorithm for blind deconvolution applied to the improvement of radiographic images, AIP Conf. Proc. 1529 (2013) 95. http: //dx.doi.org/10.1063/1.4804094.
Acknowledgments This work was partially supported by Brazilian agencies Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). References [1] W.H. Richardson, J. Opt. Soc. Amer. 62 (1) (1972) 55–59. [2] B.L. Lucy, Astron. J. 79 (6) (1974) 745–754. [3] H. Berger, Neutron Radiography. Methods, Capabilities, and Applications, Elsevier Elsevier Publishing Co., Amsterdam, 1965. [4] J.D. Domanus, Practical Neutron Nadiography, Kulwer Academic Publishers, Dordrecht, 1992. [5] J.P. Barton, Materials Evaluation (1967) 45A. [6] ASTM E803-81, Standard method for determining the L/D ratio of neutron radiographic beams.
181