=Paper=
{{Paper
|id=Vol-2485/paper5
|storemode=property
|title=RANSAC ART Tomography
|pdfUrl=https://ceur-ws.org/Vol-2485/paper5.pdf
|volume=Vol-2485
|authors=Vladimir Afanasyev,Alexey Voloboy
}}
==RANSAC ART Tomography==
https://ceur-ws.org/Vol-2485/paper5.pdf
RANSAC ART Tomography
V.V. Afanasyev1, A.G. Voloboy2
vvafanasjev@mail.ru|voloboy@gin.keldysh.ru
1
CMC MSU, Moscow, Russia;
2
Keldysh Institute of Applied Math RAS, Moscow, Russia
This paper describes using of per-voxel RANSAC approach in ART tomography. The method works as an addition to any ART and
does not depend on its internal details. Firstly, the histograms of voxel map corrections are constructed in each voxel during usual pass
of ART. Then, they are used to refine the absorption map. It allows to improve resulting voxel absorption map, reducing ghost effects
caused by input data errors and inconsistency. This method was demonstrated with optical tomography algorithm as it has certain
difficulties with input data. The proposed algorithm was implemented to run on GPU.
Keywords: Tomography, Optical Tomography, ART, RANSAC, GPU.
typical errors that lead to significant errors in final result. There
1. Introduction may be:
1. Inaccurate object models.
1.1 Tomography 2. Non-transparent regions.
Computed tomography (CT) is a class of problems of object 3. Errors in lighting model.
internal structure reconstruction using a set of its projections. 4. Errors in camera calibration.
Computed tomography methods are used in X-ray and optical 5. Noise.
tomography, depending on radiation wavelength. [1] All these types of input data inconsistency lead to same
The tomography scanner typically consists of emitter, effects on computational algorithms so we wonβt consider them
detector and a place for observable object between them. Emitter separately.
irradiates rays with known intensity, they are absorbed inside the Optical tomography is the hardest type of tomography
object and the detector registers residual ray intensity. During the problems because it has to deal with light reflection and
object rotation a set of projections (Fig. 1) is created and using refraction on all optical borders. The methods used in optical
them the internal object structure is reconstructed. tomography are applicable to X-ray tomography also.
Our work is related to optical scanning of diamonds so this
paper will focus on optical tomography. The test sample shown
here is a raw diamond molded inside a cube of immersion glass
that has index of refraction close to diamond. This helps to
eliminate scattering and refractions on rough diamond surface.
Fig. 2 shows an example of image taken in optical
tomography scanner, a cube of immersion glass with diamond
molded inside it.
Fig. 1. Parallel projection of function.
The radiation absorption inside the material obeys Beerβs
law:
π
πΌ = πΌ0 π β β«π π(π₯)ππ₯ ,
where πΌ is residual intensity received by detector, πΌ0 is initial
emitter intensity irradiated in this direction, π(π₯) is distribution
of absorption index along the ray.
The transform of volumetric distribution function into a set
of flat projections along parallel lines for all angles of rotation is
called Radon transform:
β Fig. 2. Example of optical tomography input image
π
(π, π) = β«ββ π(π + πβπ‘)ππ‘,
where π is starting point and πβ is direction vector of single
1.2 Algebraic Reconstruction Technique
line.
A theorem [3] says that under certain restrictions the inverse (ART)
Radon transform exists. One of the mostly used CT methods is ART first described
However, in most real applications it can only be computed in [2]. It is based on sequential correction of resulting function
numerically. Input data are discrete images taken around the stored in voxel map using its projections one-by-one. The input
object. The target absorption distribution function is computed in data images consist of separate pixels that correspond to rays
cubic voxel map inside the object. Also, the inverse Radon going through volume of interest. Pixel value is projection of all
transform is very unstable in case of input data errors but the absorption along this ray. Every projection value which is an
algorithm used for real computation should be stable. integral of initial function along the corresponding ray affects the
In this paper, optical tomography is used to demonstrate resulting function along the same ray according to some law. Fig.
proposed method because its input data usually have several 3 shows an example of possible absorption correction along ray
that increases its integral. Here we just add constant value to
Copyright Β© 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
every voxel along the ray. If measured value of absorption in
some pixel is π΄0 , current value of voxel map integral along the
ray is π΄ and ray length inside voxel map is π then we can add the
π΄ βπ΄
value 0 to all voxels along the ray and make the integral match
π
observed value exactly. However, such approach is impractical
as it leads to algorithm instability. Instead of this, a value less
than that can be added and observable integral along the ray gets
closer to the real measured value of this integral every time.
Then, the iterative process stops on some trigger. [1]
Fig. 6. ART iterations in single area
In this series the algorithm uses data from one complete
rotation around the object (400 directions per 360Β°). In each
image 1 or 2 different ray directions modify voxel map
depending on refraction on glass cube sides. On the directions
40-100 we can see formation of horizontal beam going through
inclusion because information from one side is mainly used.
Fig. 3. Absorption correction along the ray Then, directions 120-200 bring information from perpendicular
side and decrease absorption on the sides of the beam, amplifying
center area. But the beam sides are not removed completely.
Using information from other 2 sides does not remove these
ghost beams also. Due to the imperfection of input data described
above, more ART iterations only lead to oscillations of such
beams and they never disappear. However, photos from most
directions show that these regions containing ghost beams should
be completely transparent.
Thatβs why we propose using RANSAC in each voxel to
eliminate these negative effects.
2. Proposed method: RANSAC + ART
Our method is an add-on to any ART. It does not depend on
tomography type (X-ray or optical) and details of ART
implementation. It uses only the fact that ART iteratively refines
absorption map in each voxel and the assumption that there is no
bias in input data.
In this method, a histogram is assigned to every voxel in
addition to absorption value used in ART. Each bin of histogram
Fig. 4. Voxel cube (green) and its single horizontal cut (red) corresponds to an interval of possible absorption values and
stores number of absorption corrections made that get inside
binβs interval. An example of real voxel contents is shown on fig.
7:
Fig. 5. Horizontal layer of computed voxel absorption map
When applying corrections from different directions to voxel
map, absorption coefficient distribution is approximated. Fig. 6 Fig. 7. Absorption value and histogram in one voxel
shows ART process in small area of 32Γ32 voxels size (red
square on fig. 5) inside single voxel layer from fig. 4. It contains Here we can see that after a number of iterations most
a non-transparent inclusion in center. Brighter color means more corrections got into interval (0.7; 0.9] while the last iterations
absorption coefficient, black is 0 absorption.
gave value 0.62. It should be corrected so that it corresponds to This series shows how absorption in ghost beams around the
histogram maximum bin. inclusion is reduced to relatively low values compared to initial
The proposed algorithm works in 3 phases: ART result.
1. Pre-building voxel map with ART.
2. Usual ART algorithm with simultaneous histogram 3. Map segmentation
map construction. During this phase newly computed
If the final purpose is segmentation of absorption map, it can
absorption values add 1 to corresponding histogram
bin in each voxel. be performed in 3 ways:
1. Usual segmentation of resulting absorption map using
3. Absorption map refinement using histogram. Here
its values.
only corrections towards maximum bin are allowed in
2. Histogram map-assisted segmentation using both maps
each voxel. Corrections that make final absorption
farther from maximum bin are ignored. Corrections after absorption map refinement.
3. Segmentation using only histogram map without
inside the maximum bin can occur any number of
absorption refinement. This can significantly reduce
times, e.g., absorption can go up as far as it can in the
computation time.
last bin to represent opaque regions.
The result histogram map (after step 2) is shown on fig. 8:
4. GPU implementation
The proposed algorithm was implemented on GPU using
Nvidia Optix library [4]. Due to limited GPU memory, the
following parameters of voxel map were used:
Resolution 512Γ512Γ512
Precision Float32
Histogram bins per voxel 16
Histogram type Uint16
The algorithm consumes 5 GB of memory (10 times more
than ART) and runs on GTX Titan for 1 hour per 1200 iterations
which is about 3 times longer than ART implementation on the
same hardware.
5. Results
Fig. 8. Histogram and its binsβ intervals in voxel layer
Values of histogram are represented with RGB channels by
3 in one picture. Values 1, 2, 3 go to R, G, B of first picture; 4,
5, 6 to second picture and so on. 16 histogram bins were used in
this example. All values in this example are clamped at 255 to fit
into the image.
Fig. 9 shows magnified histogram for the same area as in Fig.
6.
Fig. 9. Magnified histogram region
Here we can see maximum in the highest bin where
πππ ππππ‘πππ β (3, β) in the center region. In this area inclusion
is located. In other bins of the histogram we can see the Fig. 11. Comparison of results
surrounding ghost beams, but unlike pure ART here they get into Fig. 11 shows the comparison of ART and proposed method
rather low bins of the histogram and their absorption will be results on the same real input data. JET color scheme is used to
significantly reduced during step 3 of the algorithm. represent details better. In the region a-b we can see reduction in
Fig. 10 shows the refinement process after forming of ghost beams around inclusions. In the region c-d phantom
histogram map. absorption computed by ART because of cube corner model
mismatch, is gathered exactly in this corner by the new algorithm
and surrounding area is almost clear as it should be, as it is in
object photos (fig. 4).
6. Conclusion
The proposed method can be used with any ART algorithm
for X-ray or optical tomography. It does not depend on such
details as reflections and refractions on optical surfaces, it is
improvement of ART itself and can be used in its other
applications. It does not require any additional restrictions on
input data compared to ART.
Fig. 10. Absorption map refinement The main advantage is better stability to input data errors
than original ART algorithms. In ART, if there is any outlier in
input data, it can damage the whole result and our method allows
to ignore such false data.
Also, small deviations in input data make ART to oscillate
around some solution so there is no sense in making high number
of iterations in order to achieve better quality. Using our method
allows to remove these oscillations.
This method works well applied to optical tomography real
data and improves results of ART algorithms, however,
theoretical question of introducing bias into result remains a
subject of further research.
The algorithm has also several disadvantages.
First, the proposed algorithm still has lack in quality, it does
not remove the ghost beams completely. However, it has not
been fine-tuned enough and has a potential.
A lot of excess memory is used to store histograms. Not
every GPU can run such algorithm with reasonable voxel map
resolution. But the amount of used memory can be adjusted by
selecting number of bins and it is not necessary to set it as high
as in this paper. Also, the memory problem can be decreased by
dumping parts of voxel map to CPU memory. Finally, in several
applications octrees and similar structures can be used to store
sparse voxel maps and reduce used memory size.
The computation time is also higher (3 times in our
implementation) due to construction and traversal of histogram
voxel map. This problem can be solved by using more than one
GPU. The algorithm can be paralleled between GPUs by splitting
voxel map and all computations along vertical axis.
7. References
[1] Afanasyev V., Ignatenko A., Voloboy A. Simultaneous
Absorption and Environment Light Reconstruction in
Optical Tomography Problem, WSCG 2015
[2] Gordon R., Bender R., Herman GT. Algebraic
reconstruction techniques (ART) for three-dimensional
electron microscopy and x-ray photography. Journal of
Theoretical Biology.29. No.3, pp. 471β81. 1970..
[3] Helgason S. Radon Transform Second Edition. Cambridge,
pp.15-20, 1999.
[4] Parker S., Bigler J., Dietrich A. OptiX: a general purpose
ray tracing engine. ACM SIGGRAPH 2010 papers, Article
No. 66. 2010.