This paper presents a multiscale feature extraction technique for 3D range images. Optimized and improved 3D Gaussian filtering and saliency map computation have been conceived jointly to the exploitation of connectivity relationships naturally induced by the 2D acquisition grid. The proposed algorithmic and implementation solutions guarantee good repeatability of detected features on different views and demonstrated superior computational performance compared to other known approaches

Alberto Signoroni

Francesco Bonarrigo

Matteo Carocci

Riccardo Leonardi

BONARRIGO, Francesco

SIGNORONI, Alberto

LEONARDI, Riccardo

Archivio istituzionale della ricerca - Università di Brescia

A MULTISCALE FEATURE EXTRACTION APPROACH FOR 3D RANGE IMAGESF. Bonarrigo∗, A. Signoroni, R. LeonardiUniversity of BresciaDept. of Information Engineeringvia Branze 38, I25123, Brescia, ItalyM. CarocciOpen Technologies srlvia IV Novembre 14I25010, S. Zeno Naviglio, Brescia, ItalyABSTRACTThis paper presents a multiscale feature extraction techniquefor 3D range images. Optimized and improved 3D Gaussianfiltering and saliency map computation have been conceivedjointly to the exploitation of connectivity relationships natu-rally induced by the 2D acquisition grid. The proposed al-gorithmic and implementation solutions guarantee good re-peatability of detected features on different views and demon-strated superior computational performance compared to otherknown approaches.1. INTRODUCTIONNowadays, there is a growing interest in 3D shape acquisi-tion and analysis because of their broad range of applicationfields (industrial, entertainment, security, medical, culturalheritage,...). In order to acquire the 3D shape of an object,a stereo scan system is usually employed to obtain a set of 3Drange images (RI). If no information about relative displace-ments between scanner and object is available, the differentviews will not share a common reference system. In order toreconstruct the object shape, it is thus necessary to accuratelyalign each view of the acquisition dataset through a registra-tion procedure. A way to estimate the relative displacementbetween two views consists into the identification of a set ofcorresponding points taken from each view and the computa-tion of an optimal rototranslation matrix out of them. Bruteforce approaches to solve this problem become impractical asRI size increases. In order to speed-up this procedure, a com-mon technique is to reduce the number of candidates by au-tomatically extracting a small set of feature points from bothviews and searching for correspondences between them. Theproblem then becomes how to select feature points that, withhigh repeatability rate, appear in same positions on overlap-ping areas of different views, i.e. taken from different angledirections. Feature extraction repeatability performance areimportant within a feature based model acquisition pipeline inthat it heavily influences subsequent steps toward automatedview alignment for 3D model acquisition [1]. Different so-lutions have been proposed for 3D feature extraction. In [2]∗The first author performed the work while at Open Technologies srl, Italypoints are selected based on curvature criteria. In [3], for eachpoint of the RI, the volume portion of the object inscribedinto a sphere surrounding the point is estimated and only un-common values are selected. Other solutions, like the pre-sented one, are based on multiscale analyses [4], [5], [6] andare somewhat inspired to the Lowe’s SIFT approach [7]. Al-though differing in algorithmic and implementation aspects, atypical sequence of processing steps for feature point extrac-tion is:• calculate the normals to the RI points (preliminary step)• calculate N Gaussian filtered (scaled) versions of the RI• derive N-1 saliency maps from the filtered data• identify local maxima on each saliency mapFollowing this approach, we present some improving solu-tions for feature detection on 3D range data and we measuretheir performance in terms of feature repeatability and com-putation speed. Our main contributions are 1) a new direction-constrained Gaussian-like filtering, and 2) naturally inducedneighborhood exploitation for computational speed-up.2. FEATURE EXTRACTION ON 3D RANGE IMAGESIn a RI, each acquired pixel (usually on a image sensing gridx = (w, h)) is associated to a depth or distance information(according to the modeled acquisition geometry). Some scan-ners produce true 3D localization of the point in the acqui-sition coordinate system. In our formalism a range imageis a map RI : I(Z2) → R(R3), where the domain I is aregular rectangular grid (usually corresponding to the CCDmatrix), and the co-domain R corresponds to the set of 3Dpoints belonging to the target object. Not all points in I havea valid corresponding point px ∈ R (this is due to the na-ture of the acquisition in terms of target object shape and/ormeasure range limitations) so a valid point subset Iv ⊆ Ican be defined. Moreover, although the 2D index grid is uni-formly distributed, 3D points cannot be considered uniformlydistributed over the 3D space. This is because point densityis inversely proportional to the angle existing between the ac-quisition direction, and the normal associated to each point.Considering R as a surface, its point distribution presentsholes in correspondence to occlusions, i.e. parts of the tar-get object that are not directly visible by the incident rays.Despite these facts (non-valid regions, non-uniform densityand holes) the regularity of the 2D grid is a useful structuringresource that can be used as a naturally induced connectivitylattice to be exploited in 3D point processing and to determinepractical definitions to various data processing operators (e.g.filtering, interpolation/subsampling, normal computation,...)and speed-up their computation. We exploit this connectiv-ity in all the steps of our method, and we propose originalimplementations of Gaussian filtering and saliency map cal-culations on RIs.2.1. Normal extractionOur feature extraction algorithm requires the computation ofnormal vectors at each point location px ∈ R. Given a px,its normal is computed by calculating the difference vectorsd = pj − px between px and its 8 neighbor points pj in the2D index grid. By assigning an index k, ranging from 1 to8, for each vector d, associated in a anti-clockwise order (seeFig.1), we define the normal nˆx as:nˆx =8∑k=1dˆk × dˆj8j ={k + 2 k <= 6k − 6 k > 6 (1)(a) (b)· · ·(c)Fig. 1. dˆ selection for (a) k = 1; (b) k = 2; (c) k = 82.2. 3D Gaussian filteringAn important step for our algorithm requires to extract a set ofN Gaussian filtered versions G (n) of the input range image,with n = [1, N ]. In order to do so, N 3D Gaussian filters withprogressively increasing kernel dimension σn are applied toeach point constituting the original range image, varying itsposition as follows:pfx (n) =∑pj∈ηx(2σn)pj · e−‖px−pj‖22·σ2n∑pj∈ηx(2σn)e−‖px−pj‖22·σ2n(2)where σn represents the nth kernel radius, while ηx (2σn)identifies the points within a distance 2σn from px.Fig. 2. pfx (n) projected over nˆx directionThe set of pfx (n), as calculated in (2), could be directly used(as done in [5] and [6]) as filtered data. However, we do notyet consider it as such. This is because, due to the above de-scribed characteristics of R, border points (on validity regionor hole boundaries) and points lying between two regions withdifferent normal directions would suffer from shrinking andposition biasing distortions respectively. These distortionscannot be considered negligible, since most feature points arelocated over highest curvature regions and therefore they areprone to position biasing due to point density disparity overηx (2σn). In addition, shrinking distortion would affect bor-der areas and would lead to artifacts in the saliency compu-tation, as we will see. In order to avoid such distortions, wepropose to project pfx (n) over the normal direction nˆx, as fol-lows:Gx (n) = px +〈pfx (n)− px, nˆx〉nˆx (3)This is equivalent to force pfx (n) to move along the axisassociated to nˆx, therefore reducing the distortion caused bynon-uniform points distribution over ηx (2σn), as depicted inFig.2.We define G (n), the set of points Gx (n) ,∀x ∈ RI, as the fil-tered version of the range image after having applied a Gaus-sian filter of radius σn, as defined in (3). As the kernel radiusσn increases, details which size is smaller than σn disappearfrom G (n). Once the G (n) has been calculated, its normalvectors are recomputed as indicated in section 2.1.As usually done in multiscale approaches, when kernel sizedoubles, a subsampling factor of two is applied to the image,also to reduce the computational burden. Here subsamplingis performed on I , thus in a very simple and fast way.2.3. Saliency maps calculationOnce the G (n) have been computed, N-1 saliency maps canbe derived. A saliency map is a 2D array of difference values,obtained by pairwise subtraction of G (n) at different kernelsizes. In a multiscale framework, a pairwise subtraction be-tween G (n) retains only the details comprised between thetwo scale limits (which can be roughly thought as cut-offfrequencies), in other words it highlights features which di-mension is comprised between two kernel sizes, say σn andσn+1. Thanks to the shrinking artifact avoidance guaran-teed by (3), straightforward pointwise subtraction at differentscales can be performed correctly. We then calculate saliencymaps S (n) as follows:S (n) = |G (n)−G (n+ 1)| · 〈nˆx (n) , nˆx (n+ 1)〉 (4)Unlike [4] and [5], the factor 〈nˆx (n) , nˆx (n+ 1)〉 has been in-troduced in order to give more importance to points for whichthe normal direction does not change after the filtering pro-cess, thus restricting saliency distribution to more specific ar-eas of the range image. A similar approach is employed in[6], where the correction factor takes into account nˆx (n).2.4. Maximum extractionOnce each saliency map S (n) has been determined, saliencymaximums are extracted as follows:sMax (n) = {x | Sx (n) > Sj (n) , ∀j ∈ ηx (2σn+1)} (5)In order to generate the maximum list for a saliency map, aniterative search is performed. Once the greatest valid saliencyvalue for S (n) is found, no other maximum may be selectedwithin a neighborhood region of ηx (2σn+1). This is becausethe greatest detail size that can be detected within S (n) isσn+1. Such constraints on the distance between maxima guar-antee the absence of overlap between features.Each maximum is further tested through two conditions be-fore being selected as feature point. At first, its neighbor-hood must be defined. In order to avoid selecting featureswhich descriptor would result calculated out of few neigh-bor points, maximums that possess less than 50% of the es-timated points for ηx (2σn) are discarded. The second checkis about position stability. Saliency maps may slightly varydue to changes of acquisition direction, therefore maximumsdetected over regions with similar saliency values may varytheir position significantly. In order to avoid such points, weonly retain maximums for which the saliency value is muchgreater than its neighboring points, and it is not located oversaliency ridges. An example of saliency map with superim-posed feature points at a given scale is shown in Fig.3.3. RESULTSTo assess the contribution of each step and the overall per-formance of the proposed feature extraction method, four dif-ferent implementations have been tested, depending on whichversion of filtering and saliency computation have been used.The first version, we called basic implementation, reflectswhat proposed in [5], i.e. Gaussian filtering performed as in(2) and saliency computed as pairwise difference between fil-tered versions of data. The second version only implementsthe improved version of Gaussian filter (3); the third ver-sion only implements the saliency correction (4). Al last, theFig. 3. A saliency map with superimposed feature points.Fig. 4. Datasets: (a) Angels, (b) Dolphin, (c) Teeth.fourth version implements both the proposed improvements.The comparison has been performed over three pre-alignedrange image datasets which are shown in Fig.4 (where differ-ent colors have been used for each view). All RIs have beenacquired with a commercial structured-light scanner (OpticalRev-Eng LE, produced by Open Technologies srl, Italy). Allthe range images used for our tests have a resolution of 1280× 1024, which means that Iv can have up to over 1.3 mil-lion points. Prior to the first Gaussian filtering, each rangeimage is down-sampled by a factor 2, thus reducing initial di-mensions to 640×512. In our implementation, three G (n)have been calculated, with σn = 2σn−1, while the σ1 is setto be eight times the average interdistance between points inR. Although, according to the SIFT methodology [7], theimplemented software can work with any number of scalinglevel and σ progression, we experimentally found the aboveparameters to give the best compromise between repeatabilityand computational performances for the considered datasets.Repeatability performance evaluations have been executed asfollows: from each range image of a dataset, its features areextracted, then for each pair of consecutive range images (ac-cording to a predetermined order), features in the overlappingarea are counted; subsequently, both feature sets are com-pared in order to find correspondences (feature points are con-sidered correspondent if their position differ in less than σ3 ).At last, feature repeatability is computed as the number offeatures that have a correspondence relation divided by theFeat. extraction versions Feature repeatabilityAngels Dolphin TeethBasic implementation 37% 8% 19%Gaussian projection 56% 34% 41%Saliency projection 35% 11% 19%Proposed implementation 57% 36% 41%Table 1. Repeatability performanceDataset RI Min-max Avg points Exec timename No. overlap (×103) [ms]Angels 8 20-60% 1040 4823Dolphin 20 20-90% 406 1894Teeth 8 40-90% 408 1836Table 2. Computational performancenumber of features in the overlapping area. In Table1, featurerepeatability performances are presented for each dataset andfor each feature extraction version. Reported data are aver-aged over all scales and all couples of overlapping RIs. Ourevaluation shows that the modified Gaussian filter improvesrepeatability performances by more than 20 percentage pointswith respect to the basic implementation of the Gaussian fil-ter (used in the first and third version). Such improvementis due to the fact that the proposed filter does not introduceany shrinking effect usually associated to the Gaussian filter-ing, thus improving feature localization. The saliency correc-tion does not yield as much, however a slight performanceimprovement is achieved through its introduction. It can beargued that such repeatability rates may not be considered asextraordinary. As a matter of fact, to our knowledge, no studyabout feature repeatability have been published. In [4], somefeature repeatability results against range data subsamplingand noise addition are presented, but they are not useful herefor comparison. We deem our performances to be satisfactorybecause of two reasons. The first one is related to the imagesets employed, since they present a heterogeneous grade ofregularity, different feature types and dimensions, and usuallybig differences (borders, holes, point density) when acquiredwith different acquisition directions. Secondly, we deem suchrepeatability rates as entirely adequate for a robust subsequentsearch of correspondences, if suitable feature descriptors orsignatures can be calculated and used to this end.Computational performances are presented in Table 2, wherewe also describe datasets in terms of number of RI views perdataset, min-max overlap area between views pairs, and av-erage number of points (i.e. |Iv|) per RI. Algorithms havebeen implemented in C++, and executed on a PC with INTELCORE2 CPU, 2.80 GHz each, and 4 GB of RAM. As thedifferences in computational time for the above four versionsare negligible, we only give the processing time required bythe proposed implementation, which is around 4.5 secondsfor each million of points. This is quite an improvement withrespect to [5] (around 5 minutes per million points), [4] (50minutes per million points), [3] (around 5 minutes per millionpoints) and [6] (around 7 minutes per million points), evenconsidering the upgraded computational resources. Such per-formance boost is mainly due to the exploitation of neigh-borhood relationships defined on I . Secondly, the use of afast subsampling, and some coding efficiency solutions onfrequent and time-consuming operations make the rest (e.g.look-up tables for square root computing).4. CONCLUSIONIn this paper a multiscale method for 3D feature extraction onrange images has been presented. This method, that can alsobe considered as a 3D extension of a SIFT-based approach,demonstrated to yield good robustness in terms of feature re-peatability, and outstanding computational performance. Thework is mainly characterized by the introduction of a modi-fied version of Gaussian filtering of 3D range images and aweighted saliency map computation. They both contributeto improve feature repeatability. In addition, the direct ex-ploitation of the 2D domain neighborhood helps the definitionof the filter action and is determinant for the computationalspeed. Experimental performance confirm that our method issuitable to be used for automated multiple view registrationand 3D model acquisition, where the degree of accuracy andobtainable automation strongly depends from the goodnessof the feature extraction phases. The proposed solutions canbe also useful for other 3D feature-based applications, suchas content analysis and description of 3D objects, 3D objecttracking, etc., where features robustness and computationalspeed are crucial aspects.5. REFERENCES[1] F. Bernardini and H. Rushmeier, “The 3D model acquisitionpipeline,” Computer Graphics Forum, vol. 21, no. 2, pp. 149–172, 2002.[2] S. M. Yamany and A. A. Farag, “Surfacing signatures: An orien-tation independent free-form surface representation scheme forthe purpose of objects registration and matching,” IEEE Trans.on Pattern Analysis and Machine Intelligence, vol. 24, no. 8, pp.1105–1120, 2002.[3] N. Gelfand, N. J. Mitra, L. J. Guibas, and H. Pottmann, “Robustglobal registration,” in Symp. Geom. Processing, 2005.[4] X. Li and I. Guskov, “Multi-scale features for approximatealignment of point-based surfaces,” in Symp. Geom. Process-ing, 2005.[5] C. H. Lee, A. Varshney, and D. W. Jacobs, “Mesh saliency,” inSIGGRAPH ’05, 2005, pp. 659–666.[6] U. Castellani, M. Cristani, S. Fantoni, and V. Murino, “Sparsepoints matching by combining 3D mesh saliency with statisticaldescriptors,” Computer Graphics Forum, vol. 27, no. 2, pp. 643–652, 2008.[7] D.G. Lowe, “Distinctive image features form scale-invariantkeypoints,” International Journal of Computer Vision, vol. 60,no. 2, pp. 91–110, 2004.

A multiscale feature extraction approach for 3D range images

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Università di Brescia