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.