Using Curvilinear Features in Focus for Registering a Single Image to a 3D Object
Hatem A. Rashwan Sylvie Chambon Pierre Gurdjos Géraldine Morin Vincent Charvillat










Using Curvilinear Features in Focus for Registering a Single Image to a 3D Object
    Hatem A. Rashwan,
        Sylvie Chambon,
        Pierre Gurdjos, Géraldine Morin and Vincent Charvillat
University of Toulouse, IRIT
    December 30, 2023
=========================================================================================================================================






In the context of 2D/3D registration, this paper introduces an approach that allows to match features detected in two different modalities: photographs and 3D models, 
by using a common 2D reprensentation. 
More precisely, 2D images are matched with a set of depth images, representing the 3D model. 
After introducing the concept of curvilinear saliency, related to curvature estimation, we propose a new ridge and valley detector for depth images rendered from 3D model. 
A variant of this detector is adapted to photographs, in particular by applying it in multi-scale and by combining this feature detector with the principle of focus curves. 
Finally, a registration algorithm for determining the correct viewpoint of the 3D model and thus the pose is proposed. 
It is based on using histogram of gradients features adapted to the features manipulated in 2D and in 3D, and the introduction of repeatability scores. 
The results presented highlight the quality of the features detected, in term of repeatability, and also the interest of the approach for registration and pose estimation. 









§ INTRODUCTION




Many computer vision and robotic applications are used to take 2D contents as input, 
but, recently 3D contents are simultaneously available and popular. In order to benefit from both modalities, 2D/3D matching is necessary. 
For medical imaging, registration of pre-operative 3D volume data with intra-operative 2D images becomes more and more necessary to assist physicians 
in diagnosing complicated diseases easily and quickly <cit.>. 
For robotic, the 2D/3D matching is very important for many tasks that need to determine the 3D pose of an object of interest: 3D navigation or object grasping <cit.>. 


The main goal of 2D/3D registration is to find the transformation of the 3D model that defines the pose for a query 2D image.
Thus, a typical 2D/3D registration problem consists of two mutually interlocked subproblems, point correspondence and pose estimation.


To match 2D photographs directly to 3D models or points clouds, most systems rely on detecting and describing features
on both 2D/3D data and then on matching these features <cit.>.

Recently, some approaches are based on learning by specific supervision classifier <cit.>. 

These methods produce very interesting results, however, they require huge amount of viewpoint-annotated images to learn the classifiers.


What makes difficulty to the problem of matching 3D features of an object to 2D features of one of its photographs is that the appearance of the object dramatically depends on intrinsic characteristics of the object, like texture and color/albedo, as well as extrinsic characteristics related to the acquisition, like the camera pose and the lighting conditions.




Consequently, some approaches manually define correspondences between the query image and the 3D model, such as <cit.>.
These manual methods can be robust but it can easily become hard to apply this manual selection to large image sets. Moreover, in this paper, we focus on automated approaches.


Note that some systems are able to generate a simultaneous acquisition of photographs and scanning of a 3D model but using this kind of systems induces limited applications.

Other methods solve the problem by distinguishing two subproblems: to choose the common representation of the data and, then, to find the correspondences. 
These methods transforms the initial 2D/3D registration problem to a 2D/2D matching problem by rendering multiple 2D images of 3D models from different viewpoints, 
such as <cit.>.

Consequently, the first task of 2D/3D registration is to find an appropriate representation of 3D models in which reliable features can be extracted in 2D and 3D data. 
In <cit.>, synthetic images of the 3D model are rendered, while depth images are rendered in <cit.>. 
More recently, <cit.> proposes average shading gradients. This rendering technique for a 3D model averages the gradient normals over all lighting directions 
to cope with the unknown lighting of the query image. 
The advantage of representing the 3D model by a set of depth images is that it can express the model shape independently to color and texture information.
Therefore, representing the 3D model by a set of depth images is the best option for this work, see Fig. <ref>. 
In this case, features extracted from depth images are only related to shape information.




















Then, the second difficulty of 2D/3D registration consists in proposing how to match entities between the two modalities in this common representation. 
It can be partial <cit.> or dense matching, based on local or global characteristics <cit.>. 
In <cit.>,  silhouettes extracted from synthetic images are matched to ones extracted from the color images. 
However, this method is not able to take into account most of the occluding contours that are useful for accurate pose estimation. 
In turn, in <cit.>, image gradients are matched with their 3D representation. 
Since image gradients are still affected by image textures and background, this technique can fail to estimate the correct correspondences. 
A key requirement on these features, as in classic 2D matching between real images, is to be computed with a high degree of repeatability. 
In our case, similar to the definition in <cit.>, 
the repeatability of a feature is defined as the frequency with which one detected in the depth image is found 
within ϵ pixels around the same location in the corresponding intensity image (if it is supposed that the features are not moving or are following a small deplacement).

Then, since we suppose that an individual photograph of an object of interest is acquired, in a textured environment, we will focus on comparing pre-processed features of color images
with features of a set of rendered images of 3D models, more precisely, a set of depth images, see Fig. <ref>.

More precisely, the 3D object will be given by a set of 3D depth surfaces, which describe how the original object surface is shortened by a perspective viewing and 
the image is given by the 3D intensity surface. Since the depth and the intensity surfaces have a different order of representation, the two surfaces can not be directly matched.

Thus, bringing both rendered depth images and photographs into a common representation, such as gradient and edge representation, allows to establish a robust sparse 2D-to-3D matching <cit.>.

We propose to extract gradient-based features corresponding to object's shapes in both depth and intensity images regardless of illumination and texture changes.
In other words, as 2D photographs (intensity images) are affected by background, textures and lighting changes,
we take into account these difficulties by reducing the influence of non-redundant information (i.e., color and texture) on features extracted from photographs.


It means that we extract features in depth images that highlight geometric characteristics of an object.
For photographs, we need to refine detected features by selecting salient points acquired by a camera in focus.
These points are a function of the degree of blur (blurriness) in an image.
Thus, the detected points are analyzed based on measuring the blur amount of every feature point.
Finally, what we call focus points should be able to detect the approximate shape and to discard the other components such as textures.



To summary, the contributions of this paper, as shown in Fig. <ref>, are:

 
  * A ridge and valley detector, for depth images rendered from 3D model. We name it curvilinear saliency (CS) as it is related to the curvature estimation
 (a function of the eigenvalues derived of the Hessian matrix).
 This representation directly relates to the discontinuities of the object's geometry, and, by nature, the extracted features should be robust to texture and light changes.

  * A variant of this detector adapted to photographs.   
This curvilinear saliency detector is applied in multi-scale by searching over all scales and all image locations
in order to identify scale-invariant interest points.
To reduce the influence of structures due to texture and background regions, we introduce the extraction of focus curvilinear saliency features.
It corresponds to ridges that are not affected by blur.


  * A registration algorithm for determining the correct viewpoint of the 3D model and thus the pose.  
This method is based on using histogram of gradients, HOG, features <cit.>, adapted to the features manipulated in 2D and in 3D, and the introduction of 
repeatability scores.
More precisely, the HOG descriptor is computed on both depth images (i.e., curvilinear features extracted with curvilinear saliency detection) and photographs (i.e., curvilinear features in focus extracted with multi-scales curvilinear saliency detection) and it combines the curvilinear saliency value with the orientation of the curvature. The repeatability score measures the set of repeatable points detected both in a photograph and in the rendered depth images.


After presenting the related work and reminder on differential geometry, sections <ref> and <ref>, 
we introduce the 3D model reprensentation, section <ref>, and, then, the image representation used for 2D/3D matching, section <ref>. 
Then, we describe how we try to be robust to background and texture by using the same principle used in the detection of focus curves, section <ref>. 
We illustrate how this new global approach for 2D/3D matching allows to obtain more repeatable features, compared to state of the art, section <ref>. 
Finally, we explain how we obtained 2D/3D registration results, section <ref>, 
and pose estimation, section <ref>, by highlighting the interest of the proposed approach in these applications, before conclusion, section <ref>. 



§ RELATED WORK

As mentioned earlier, a typical 2D/3D registration problem consists of two subproblems: feature correspondence and pose estimation (i.e., alignment). Thus, the related work is divided into three parts related to these subproblems: 1) detect features in 2D photography 2) detect features in a 3D model and finally 3) match 2D features to 3D features to estimate the 3D pose.




 §.§ Classic 2D feature detection





In this section, we try to identify if, in the literature, it exists a 2D classical detector interesting in order to obtain points comparable as points detected, 
with the same principle or same tool, in 3D. 
We suppose that, for this purpose, it is necessary to detect features that are related to points of interest on the structure of the object and not on the texture or the 
light changes on the object. 
In 2D, edge detection <cit.> based on the first-order derivative information is the initial technique. 
It can detect any kind of edges, even low contrasted edges that are not due to the structure but more due to texture. 

The second technique is to detect the interest points <cit.> that refer to point-like features in an image by analyzing the eigenvalues of the structure tensor 
(i.e., the first-order derivative) at each point – the two eigenvalues have to be maximal to highlight a point of interest.
Again, this technique does not take into account difficulties due to textures, light changes or scale changes. 

In another way, blob detection <cit.> provides a complementary description of image structures in terms of regions, as opposed to point-like algorithms.
These methods are based on the Hessian matrix (i.e., the second-order derivative) because of its independence to zero- and first-order changes and its good performance in computation time and accuracy.



More recently,  multi-scale approaches have been introduced, like a generalization of Harris or Laplacian detectors <cit.> or
the well known approach of SIFT, Scale Invariant Features Transform <cit.>.
In <cit.>, SURF, Speeded Up Robust Features, a detector also based on Hessian matrix analysis, is introduced to be faster than SIFT and other multi-scale techniques by using approximation of Laplacian and algorithmic tricks.
All these techniques are robust to light changes, rotations and translations. It makes the features detection invariant to view-point changes.

However, they totally rely on texture and/or intensity changes to find the features.

Curvature detection is one of the most important techniques of second-order derivative-based approaches used for extracting the structure properties. 
Recently, <cit.> has proposed a detector based on curvature κ expressed as the change of the image gradient along the tangent to obtain a scalar q 
approximating κ. 
In addition, <cit.> presented PCBR, Principal Curvature-Based Regions, 
detector that uses the maximum or minimum eigenvalue of the Hessian matrix to find the principal curvature in a multi-scale space. 
As mentioned in <cit.>, the maximum eigenvalue yields a high value only for the dark side of edges, 
i.e. the minimum eigenvalue detects light lines on a dark background. 
By definition, it is a restrictive way to select features and it does not guarantee to select the maximal number of reliable features.

In conclusion of this review, curvature features have several advantages over more traditional intensity-based features <cit.>, 
especially with extracting local structure of the points of interest. 
In addition, particularly, curvature features are invariant to viewpoint changes and to transformations that do not change the shape of the surface.





 §.§ Classic 3D feature detection

In this subsection, our goal is to find if, in the literature, some 3D detector can be directly adapted to 2D, in order to obtain comparable points of interest.

Feature extraction of 3D models/scenes can be classified into point-based and image-based approaches.
Most of point-based methods are based on using SIFT in 3D by proposing an adaptation of the initial SIFT <cit.>.
These approaches are interesting but they are not dedicated to 2D/3D registration and, so, they do not consider to detect similar features in 2D and in 3D.



Other methods proposed to find curves that have special properties in terms of differential geometry of the surface.
For example, in <cit.>, curves are supposed to be located at paraboloic lines which occur at points of extremal curvature.
These curves capture important object properties closed to the object surface, but they do not vanish along the surface when the viewpoint changes.









In image-based approaches, the 3D model is first rendered to form images or geometric buffers.
Then image processing methods are applied such as edge detection <cit.>.
In <cit.>, 2D SIFT is applied on images of a rendered 3D mesh because this multi-scale representation extracts features 
that are supposed to be related to local extrema of the surface variation.
The apparent ridges (AR), which are a set of curves whose points are local maxima on a surface, are introduced in <cit.>.
In this paper, a view-dependent curvature corresponds to the variation of the surface normal with respect to a viewing screen plane.
However, apparent ridges often produce false edges that are not related to occluding contours, which are important for pose estimation of most of objects.
Mesh saliency measures the region importance of 3D models using Gaussian-weighted mean curvatures in multi-scales <cit.>.
However, mesh saliency is based on mesh images that are also affected by lighting conditions.
Average Shading Gradients, ASG, was proposed in <cit.>. 
This rendering technique is based on averaging gradients over all lighting directions to cope with the unknown lighting conditions. 




In conclusion, only <cit.> can be used in a 2D/3D registration context and we will describe, in section <ref>, and compare our work with this method, 
in section <ref>. 




 §.§ 2D/3D registration and pose estimation










In the computer vision literature, the problem of automatically aligning 2D photographs with an existing 3D model of the scene has been investigated in depth 
over the past fifteen years. 
In the general case, the proposed solution will be an image-to-model registration to estimate the 3D pose of the object. 
The 2D-to-3D registration problem is approached in the literature through indirect and direct methods <cit.>.

For indirect registration, these methods are performed either by 3D-to-3D registration or by finding some appropriate registration parameters, 
such as the standard Iterative Closest Point, ICP algorithm <cit.>. 


This kind of techniques are more global and do not really provide points to points correspondences. 

For direct registration methods, in <cit.>, correspondences are obtained by matching SIFT feature descriptors 
between SIFT points extracted in the images and in the 3D models. 
However, establishing reliable correspondences may be difficult due to the fact that the set of points in 2D and in 3D are not always similar, 
in particular because of the variability of the illumination conditions during the 2D and 3D acquisitions. 
In the same context, in <cit.>, the authors assume that the object in the input image has no or poor internal texture.
Methods relying on higher level features, such as lines <cit.>, planes <cit.> and building bounding boxes <cit.>,
are generally suitable for Manhattan World scenes and hence applicable only in such environments.
Skyline-based methods <cit.> as well as methods relying on a predefined 3D model <cit.> are, likewise, of limited applicability.
Recently, the histogram of gradients, HOG, detector <cit.> or a fast version of HOG <cit.> have been also used to extract the features from rendering views and real images.
All of these approaches give interesting results, however, they do not evaluate the repeatability between the set of points detected in an intensity image and those detected in an image rendered from the 3D model.
Finally, in <cit.>, 3D corner points are detected using the 3D Harris detector and the rendering average shading gradients images on each point.
For a query image, similarly, corner points are detected in multi-scale. 
Then, the gradients computed for patches around each point is matched with the database containing average shading gradient images using HOG descriptor.
This method still relies on extracting gradients of photographs affected by textures and background and it can give erroneous correspondences.
Consequently, they propose a refine stage based on RANdom SAmple Consensus, RANSAC <cit.> to improve the final pose estimation.



























In this paper, structural cues (e.g., curvilinear shapes) based on curvilinear saliency are extracted instead of only considering silhouettes,
since they are more robust to intensity, color, and pose variations.
In fact, they have the advantage to both represent outer and inner (self-occluding) contours that also characterize the object and that are useful for estimating the pose.
In order to merge in the same descriptor curvilinear saliency values and curvature orientation, the HOG descriptor, 
widely used in the literature and that properly describes the object shape, is employed.
Finally, HOG features and repeatability scores are used to match the query image with a set of depth images rendered from a 3D model.


In the rest of the paper, after the reminder on differential geometry, we will describe the 3D model representation and the image representation that are introduced in order to compare and to match the 3D and the 2D data. 
In particular, we will explain how these representations allow to be robust to background details and to texture before illustrating the interest of the proposition with 
experimental results on feature detection, registration between 2D images and 3D models and pose estimation. 



§ REMINDER ON DIFFERENTIAL GEOMETRY




 §.§ Notations





In the sequel, these notations are used: 
 
  
  * ∇_ f: the gradient vector of a scalar-valued function f. 
  
  * 𝐅_x:  the partial first-order derivative ∂𝐅/∂ x of a vector-valued
function 𝐅 w.r.t. variable x. 

  * Similarly,  𝐅_xy: the partial second-order derivatives  ∂^2𝐅/∂
x∂ y  of  𝐅 w.r.t. variables x and y.






Moreover, we assume a calibrated
perspective camera where the image
point coordinates are given  with respect to the normalized image
 frame i.e., as 𝐱 = (x,y) obtained from the equation
 [x,y,1]=𝙺^-1[u,v,1], where (u,v) are
 pixel coordinates and 𝙺 is the usual calibration upper triangular matrix <cit.>.




 §.§ Definitions: Differential of a Map – Tangent Plane – Gauss Map


Let  ℱ⊂ℝ^3 be a
regular surface[See the definition of a regular surface in R^n in <cit.>. ]

whose parameterization is given by the differentiable map
𝐅:U⊂ℝ^2→ℱ with 

    𝐅(x,y)=[X(x,y),Y(x,y),Z(x,y)]^⊤

To each
𝐱=[x,y]^⊤∈ U is associated a  map   d𝐅_𝐱:ℝ^2→ℝ^3,
called the differential of 𝐅 at 𝐱 and  defined as follows <cit.>.
Let 𝐯∈ℝ^2   be a vector and let 
α:(-ϵ,ϵ)→ U
  be a
differentiable  curve satisfying α(0) = 𝐱 and α'(0) = 𝐯.
By the chain rule, the  curve 
𝐅∘α 
in ℝ^3 is  also differentiable.  We then define
    d𝐅_𝐱( 𝐯)=(𝐅∘α)'(0)







It provides a linear (i.e., first-order) approximation 
to 𝐅(𝐱+𝐯)
when the increment 𝐯
is small enough.
This is illustrated in Fig. <ref>(a).
The vector subspace d𝐅_𝐱(U)⊂ℝ^3

 has dimension 2 and is a plane consisting of all tangent vectors of 𝐅
at 𝐏=𝐅(𝐱).
It is called the tangent plane
 of ℱ at 𝐏 and  denoted by
 T_𝐏(ℱ).

It can be proved <cit.> that the above
definition does not depend on the
choice of  α.
Furthermore, the fact that 
(𝐅∘α)'(0)=𝐅'(α(0))α'(0)
entails that d𝐅_𝐱( 𝐯) is   linear
in 𝐯. In particular, in the  canonical bases of ℝ^2 and
ℝ^3, we have

    d𝐅_𝐱(𝐯)
    =𝙹_𝐅(𝐱)𝐯

involving the 3× 2 Jacobian
 matrix 
of 𝐅 at 𝐱   

    𝙹_𝐅(𝐱)=
    [ 𝐅_x(𝐱) 𝐅_y(𝐱) ]








 with  
𝐅_x and 𝐅_y as columns. This also shows that the vector subspace d𝐅_𝐱(U)

 has indeed dimension 2.

Let  𝐏=𝐅(𝐱) be a point of ℱ. Let 𝐍:ℱ⊂ℝ^3→Σ⊂ℝ^3 be
        the differentiable map that assigns  to 
𝐏 the
 coordinate vector  𝐍(𝐏) on the unit sphere Σ representing
   the unit normal
of ℱ at 𝐏 and computed as

    𝐍(𝐏)=𝐅_x(𝐱)×𝐅_y(𝐱)/‖𝐅_x(𝐱)×𝐅_y(𝐱)‖ with 𝐱=𝐅^-1(𝐏)







    This map is called the  Gauss map of  ℱ.

The Gauss map is a mapping between the two surfaces  
ℱ and  Σ and the definition of differential 
is extended to that case.
The differential
of the Gauss map of  ℱ at point 𝐏 is the  map   
d𝐍_𝐏:T_𝐏(ℱ)⊂ℝ^3→ℝ^3
  defined as follows.
Let 𝐕∈ T_𝐏(ℱ)   be a vector on the tangent plane
 of ℱ at 𝐏 and let 
β :(-ϵ,ϵ)→ℱ
  be a
differentiable  curve  on the surface  ℱ satisfying β(0) = 𝐏 and β'(0)
= 𝐕.
By the chain rule, the  curve 
𝐍∘β 
in ℝ^3 is  also differentiable;  we then define
    d𝐍_𝐏( 𝐕)=(𝐍∘β)'(0)

It expresses how 𝐍  behaves
 — how  ℱ curves— in the vinicity of 𝐏.
This is illustrated in Fig. <ref>(b).
Again, it can be proved <cit.> that the above
definition of  d𝐍_𝐏
does not depend on the
choice of one possible curve   β.
 



Similar to the differential of a map seen above, the fact that 
(𝐍∘β)'(0)=𝐍'(β(0))β'(0)
entails that d𝐍_𝐏( 𝐕) is   linear
in 𝐕=β'(0). The vector subspace d𝐍_𝐏(T_𝐏(ℱ))⊂ℝ^3

 has dimension 2: it is the plane  
 T_𝐍(𝐏)(Σ) consisting of all tangent vectors  to the unit
 sphere at point 𝐍(𝐏) called the tangent plane of Σ at 𝐍(𝐏).
Therefore the domain of values of d𝐍_𝐏 is d𝐍_𝐏(T_𝐏(ℱ))=T_𝐍(𝐏)(Σ). Actually, it can readily seen that T_𝐍(𝐏)(Σ) and T_𝐏(ℱ) are parallel planes so the differential   of 




𝐍 is usually defined as  d𝐍_𝐏:T_𝐏(ℱ)→
T_𝐏(ℱ). 
 




 §.§ Curvatures. Fundamental forms of a surface

Let 𝐓∈ T_𝐏(ℱ)   be a unit vector representing a direction on the tangent
plane
 of ℱ at  𝐏.
Let
 𝒞  be the  curve obtained by slicing
ℱ
 with the  normal section of 
ℱ
   at 𝐏 along 
𝐓   i.e.,[This holds for any plane through 𝐏
 parallel
to    𝐓.
 This is due to
 Meusnier's theorem <cit.>  “All curves lying on a surface
𝒮
 and having at a given point 𝐏∈𝒮 the
 same tangent line have at this point the
 same normal curvatures. ”]  the plane through 𝐏  parallel to
both
𝐍(𝐏) and   𝐓.
 The normal curvature  of ℱ in   (unit) direction
𝐓
is  the curvature of  𝒞  at  𝐏   which can be given
by <cit.>: 

    κ_N(𝐓)
    =-𝐓·d𝐍_𝐏(𝐓)






It expresses how 𝐍  behaves
 — how  ℱ curves— in the vinicity of 𝐏.
An important remark is that the radius of curvature of  𝒞   at
 𝐏 is  equal to 
1/|κ_N|.

The principal curvatures  κ_1, κ_2 of ℱ
  at 𝐏 can
be defined as the extrema of  function (<ref>) 
with respect to directions 𝐓∈ T_𝐏(ℱ), subject to the constraint ‖𝐓‖=1.
 The corresponding
directions are called
principal directions  of ℱ at 𝐏.
It is well-known that ``the coefficients κ_1, κ_2 are  decisive
parameters that fully
describe local surface shape up to the second order
modulo a rigid movement.” <cit.>
This means that the principal curvatures are invariant to the surface parameterization.


Now consider a new 3D coordinate system with
𝐏=𝐅(𝐱) as  origin.



Any vector 𝐓 on  the tangent plane T_𝐏(ℱ) can
be written as

    𝐓=u𝐅_x(𝐱)+v𝐅_y(𝐱)=𝙹_𝐅(𝐱)[ u; v;   ]

where 𝙹_𝐅(𝐱) is the Jacobian matrix of 𝐅 defined in (<ref>) and (u,v) are  so-called local coordinates of  𝐓.
From now on, we will put uv as subscript to relate a vector to its local
coordinates, e.g.,   𝐓_uv.  



  §.§.§ First fundamental form of a surface

Given
any (u,v), the norm of any vector  𝐓_uv on  T_𝐏(ℱ)
writes





    ‖𝐓_uv‖=√([u,v]𝙸_𝐏[u,v]^⊤)

where

    𝙸_𝐏   =
    [ 𝐅_x(𝐱)·𝐅_x(𝐱) 𝐅_x(𝐱)·𝐅_y(𝐱); 𝐅_x(𝐱)·𝐅_y(𝐱) 𝐅_y(𝐱)·𝐅_y(𝐱) ]




The quadratic form on  T_𝐏(ℱ) 

    ℑ_𝐏(u,v)≐[u,v]𝙸_𝐩[u,v]^⊤

















 is called the   first fundamental form of ℱ
<cit.>.





























































































































































































































  §.§.§ Second fundamental form of a surface

 Let  𝐓_uv be a direction of  the tangent plane T_𝐏(ℱ), given in local 2D coordinates.  On the one hand, it can be shown <cit.> that the differential
of the Gauss map of  ℱ at 𝐏=𝐅(𝐱) along   𝐓_uv
writes in standard 3D coordinates

    d𝐍_𝐏(𝐓_uv)=[𝐍_x(𝐱)|𝐍_y(𝐱)][ u; v;   ]







On the other hand,

let denote by d𝐍_𝐏(u,v) 
the differential  of 




𝐍 at 𝐏=𝐅(𝐱) along   𝐓_uv expressed in local 2D coordinates
i.e., such that d𝐍_𝐏(𝐓_uv)
=𝙹_𝐅(𝐱)]d𝐍_𝐏(u,v).
Then we have

    d𝐍_𝐏(u,v)=𝙸_𝐏^-1𝙸𝙸_𝐏[ u; v;   ]

where

    𝙸𝙸_𝐏   =
    [ 𝐍_x(𝐱)·𝐅_x(𝐱) 𝐍_x(𝐱)·𝐅_y(𝐱); 𝐍_x(𝐱)·𝐅_y(𝐱) 𝐍_y(𝐱)·𝐅_y(𝐱) ]


The proof can be found in <cit.>.

















The quadratic form

    ℑℑ_𝐏(u,v)=[u,v]𝙸𝙸_𝐏[u,v]^⊤

















 is called the   second fundamental form of ℱ
<cit.>.
It directly follows from this that Eq. (<ref>)  can be expressed in local coordinates and writes

    κ_N(u,v)=[u,v]𝙸𝙸_𝐏[u,v]^⊤

For any (non unit)  vector 𝐕(u,v) on
 T_𝐏(𝒮), in order to set it as a unit direction, it is needed to divide (<ref>)
by the square of expression  (<ref>). Hence, the normal curvature
in  the unit direction 𝐕/‖𝐕‖
 is can be given by

    κ_N(u,v)=-𝐕_uv·d𝐍_𝐏(𝐕_uv)/‖𝐕_uv‖^2=ℑℑ_𝐏(u,v)/ℑ_𝐏(u,v)





 

  §.§.§ Closed-form solutions for principal curvatures

The principal curvatures  κ_1, κ_2 of 𝒮   at 𝐏 can
be defined as the extrema of  function (<ref>) with  (u,v)-coordinates
as variables.
 

Seeing (<ref>) as a generalized Rayleigh quotient, it
is known <cit.> that κ_N has an extremum
at (û,v̂) only if κ_N(û,v̂)
is a root of (^𝙸𝙸_𝐏-κ_N𝙸_𝐏)
 or, equivalently, only if κ_N(û,v̂) is an eigenvalue
of the 2× 2 matrix  𝙸_𝐏^-1𝙸𝙸_𝐏, which is not
symmetric but always
 has real eigenvalues <cit.>. 
As a result, the principal curvatures are the two eigenvalues  κ_α   (α=1,2) of the matrix 𝙸_𝐏^-1𝙸𝙸_𝐏. The principal 3D directions are  
𝐓_α=𝙹_𝐅(𝐱)𝐞_α where  𝐞_α  are  the corresponding eigenvectors. 

Now we state a proposition that we derive from  the above results, which will be used in our work.


The principal curvature κ_α   (α=1,2) at 𝐏
 associated to  the unit principal 3D direction   
𝐓_α

 is equal to the absolute magnitude
of the differential of the  Gaussian
map at this point   i.e.,

    |κ_α|
    =‖d𝐍_𝐏(𝐓_α)‖




Since the Euclidean norm is invariant to changes of Euclidean coordinates, without loss of generality, let choose a new parameterization
𝐒̃(x̃,ỹ)=[ x̃,ỹ,Z̃(x̃,ỹ) ]^⊤, for some  new height function Z̃, w.r.t. 3D orthonormal  frame whose origin is   𝐏 and  x̃ỹ-plane
coincides with
the tangent plane T_𝐏(ℱ).






More generally,  we will add the symbol  to distinguish the new representations from the old ones, except for the principal curvatures which are irrespective of  parameterizations. Let us remind that 
𝐓̃_α=𝙹_𝐅(𝐱)ẽ_α,
where ẽ_α is the associated eigenvector, and note that the new    first fundamental matrix,
    Ĩ_𝐏, is then the identity. As a result, starting from the fact that  κ_α is an eigenvalue
of the 2× 2 matrix  𝙸_𝐏̃^-1𝙸𝙸_𝐏̃, we have:

    (𝙸_𝐏̃^-1𝙸𝙸_𝐏̃)ẽ_α   =κ_αẽ_α
    ⇔dÑ_𝐏̃
    (ẽ_α)
       =κ_αẽ_α
    ⇔[𝐒̃_x̃|𝐒̃_ỹ]dÑ_𝐏̃(ẽ_α)
       =κ_α[𝐒̃_x̃|𝐒̃_ỹ]ẽ
    
    
    
    
    ⇔dÑ_𝐏(𝐓̃_α)   =κ_α𝐓̃_α
    ⇒‖dÑ_𝐏(𝐓̃_α)‖^2
       =κ_α^2






§ 3D MODEL REPRESENTATION





The work the most related to what is proposed in this paper is the Average Shading Gradient (ASG) approach,  proposed in <cit.>. 
After introducing how object surface can be represented, we highlight the differences between these two approaches. 



Object surface : 
Denote by   ℳ   the surface of some observed object associated to a parameterization  𝐌(x,y)≜[X(x,y),Y(x,y),Z(x,y)]^⊤,
  where (x,y) varies over
the restricted image domain of a given camera which is delimited by the occluding contour  of
 the object.
Under perspective projection, every visible 3D point of  ℳ
(seen from the camera viewpoint), with vector
𝐌(x,y), is assumed to be in one-to-one correspondence
with the 2D image point with vector 𝐱=[x,y]^⊤,
such that x=X(x,y)/Z(x,y) and y=Y(x,y)/Z(x,y).
As a result, we get
 
    𝐌(x,y)=Z(x,y)[x,y,1]^⊤

Let 𝐍(x,y) denotes the  Gaussian map of ℳ  which  assigns, on the
unit sphere, to each point point
𝐌(x,y) of  ℳ the unit normal
of ℳ defined by 
𝐍(x,y)=𝐍̅(x,y)/‖𝐍̅(x,y)‖ where, using (<ref>),
 𝐍̅ writes

    𝐍̅
    =
    𝐌_x×𝐌_y 
    
    =Z[ -Z_x, -Z_y, xZ_x+yZ_y+Z ]^⊤

It can be shown that the  Jacobian 3× 2 matrix of 𝐍 writes

    𝙹_𝐍=
    [ 𝐍_x 𝐍_y ]
    =( 𝙸-𝐍𝐍^⊤)𝙹_𝐍̅

where  the columns of 
𝙹_𝐍̅=
[ 𝐍̅_x 𝐍̅_y ]
have the form

    𝐍̅_⋆=[                     Z_xZ_⋆-Z_x⋆Z;                    Z_⋆Z_y-Z_⋆ yZ; xZ_x⋆Z+yZ_⋆ yZ+Z_⋆(xZ_x+yZ_y+3Z);                                  ]

⋆ standing for either x and or y.
 

























 §.§ Average Shading Gradient (ASG) Feature <cit.>

Plötz et al. assumed in <cit.>
that the image intensity function
obeys the Lambertian shading function


    I(x,y)∝max(0,-𝐍(x,y)·𝐬)


for a parallel light source 𝐬∈ℝ^3.




This means that the reflectance  describing the object material is assumed
 to be Lambertian with constant
albedo[A general shading function is
I(x,y)=ρ(𝐌(x,y))max(0,-𝐍(x,y)·𝐬)
where ρ(𝐌(x,y)) is the albedo at object  point 𝐌(x,y).
].
In addition, the background is assumed to be constant (e.g., a plane).



The authors propose as feature in the intensity image the magnitude of the gradient of the shading function.  
To register the intensity image to the 3D (untextured) model,  the idea is to generate virtual images viewing the object from different camera pose candidates. 
Nevertheless, it is clearly impossible to render any such virtual image obeying the shading function (<ref>) 
without prior information about the lighting direction and so about 𝐬.
Therefore, the authors propose to  replace the gradient magnitude feature, in the virtual images, 
by a feature corresponding to the average value of the gradient magnitude computed over all light directions, so-called average shading gradient magnitude.  
Denoting ‖∇_ I‖^2 = I_x^2+I_y^2 the magnitude of the gradient of the shading function (<ref>) 
then the magnitude of the average shading gradient is:

    ‖∇_ I‖=
    ∫_𝒮‖∇_I‖d𝐬

where the  vector 𝐬, cf. (<ref>),
varies over the unit sphere 𝒮 in ℝ^3 and d𝐬
is the volume element.


The nice contribution of Plötz et al.  is, by applying Jensen's
inequality, to derive the following closed-form bound on ‖∇_
I‖

    ‖∇_ I‖   ≤√(∫_𝒮‖∇_I‖^2 d𝐬)
       =γ√((‖𝐍_x‖^2+‖𝐍_y‖^2)
    )

with γ=√(π/3).
It is reported by the authors to behave like a very good approximation of
‖∇_ I‖. This is the elegant way
 the authors get rid of the unknown lighting direction
 𝐬 in (<ref>).








 §.§ Curvilinear Saliency Features (CS)


As already mentioned, our goal is to find a common representation between the 3D model and the 2D image in order to be able to match them. 
For that purpose, we first show how the 3D model can be represented or studied from different points of view and 
how these different viewpoints can be analyzed and compared to a 2D image. 
For that purpose, we   represent the observed 3D object  by a set of  synthetic
depth maps generated from  camera locations distributed on concentric spheres encapsulating, by sampling elevation and azimuth angles, as well
as distances from the camera to the object. 
A depth map (or  depth image)  Z(x,y) associates
to every image point (x,y) the Z-coordinate, w.r.t. the camera frame,
of
the object 3D point (<ref>) that projects onto (x,y).

Let 𝒟   denote the depth surface  that is the 3D surface whose
graph parameterization is[Note the difference with (<ref>).
]
 
    𝐃(x,y)=[ x,y,Z(x,y) ]^⊤







Which features should be extracted in the depth map? We aim at detecting depth
“discontinuities”  by searching points on 𝒟  having high principal
curvature in one direction
and low principal curvature in the orthogonal direction.  
We call Curvilinear Saliency features of a surface
  loci of such points. Basically, they correspond to the ridges
and valleys of this surface.
In this work, we use the difference of the principal
curvatures κ_1 - κ_2 to describe the ridges and
valleys.













Principal curvatures and directions: 


Consider  a point 𝐏=𝐃(x,y).

 Let 𝐍'(x,y) denote the  Gaussian map 
 of 𝒟
assigning to 𝐏 the unit normal of  𝒟 at 𝐏,
such that

    𝐍'=𝐍̅'/‖𝐍̅'‖ where 𝐍̅'
    =
    𝐃_x×𝐃_y
    =α[ -∇_Z;    1;      ]

with ∇_ Z=[Z_x,Z_y]^⊤ and α =1/√(1+‖∇_ Z‖^2).



 As the two columns of the Jacobian matrix
𝙹_𝐃 of 𝒟
are 
   𝐃_x=[1,0,Z_x]^⊤
and 𝐃_y=[0,1,Z_y]^⊤, the first fundamental form of 𝒟
can be computed as

    𝙸_𝐏
    =𝖨_3+∇_ Z∇_ Z^⊤



and the   second fundamental form of 𝒟
can be computed as


















    𝙸𝙸_𝐏
    =α𝙷_Z











where 𝙷_Z
is the Hessian matrix of Z i.e., with the second-order partial derivatives of Z w.r.t. x and y as
elements.











The principal curvatures 
of  𝒟      at 𝐏 coincide with  the eigenvalues











κ_α (α=1,2) of 𝙸_𝐏^-1𝙸𝙸_𝐏, which are always real. In the tangent plane T_𝐏(𝒟), 
the local coordinates of the  principal directions of 𝒟 at 𝐏 are given by the  eigenvectors 
𝐞_α of 𝙸_𝐩^-1𝙸𝙸_𝐩 so the 3D principal directions in 3D wrote 𝙹_𝐃𝐞_α. 
As Koenderink wrotes in <cit.>, 
“it is perhaps not superfluous to remark here that the simple (eigen-)interpretation in terms[By neglecting 𝙸_𝐏.] 
of 𝙸𝙸_𝐏=α𝙷_Z is only valid in representations where ∇_Z=0”, 
which is the condition for the point to be local extremum.

Thanks to proposition <ref>, we know that that the principal curvature κ_α   at 𝐏  associated to the principal 3D direction   
𝐓_α=𝙹_𝐃𝐞_α is equal to the absolute magnitude of the change of the  normal

    |κ_α|=‖d𝐍'_𝐏(𝐓_α) ‖

where  d𝐍'_𝐏(𝐓) denotes the  differential  of  𝐍' at 𝐏 in  direction  𝐓.
We will make us of this result for the image representation, cf. <ref>. 
Now let us explain why we propose as feature the difference κ_1-κ_2 where κ_1≥κ_2.
































Curvilinear feature:
Without losing generality, let κ_1 and κ_2 be the  principal curvatures computed as ordered eigenvalues of 𝙸_𝐩^-1𝙸𝙸_𝐩
so that κ_1≥κ_2. 
We aim at detecting  points lying on “elongated” surface parts.
In this work, we detect points at which this difference is  high: 

    CS(x,y)=κ_1(x,y)-κ_2(x,y)

We call (<ref>) the curvilinear saliency (CS) feature. Curvilinear means a feature that belongs to a curved line. 
The rest of this paragraph justifies such a choice.

Given a point 𝐏 on 𝒟, let (x̃,ỹ) be the Cartesian coordinates on the tangent plane T_𝐏(𝒟))
w.r.t. the 2D  frame whose origin is 𝐏 and the orthonormal basis is formed by the principal directions {𝐞_1,𝐞_2}. 
As a result, 𝒟 can now locally be associated to the new parameterization 𝐅(x̃,ỹ)=[ x̃,ỹ,F(x̃,ỹ) ]^⊤, 
for some  height function F. 
In that case, it can be readily seen that 𝙸_𝐏 is the identity matrix, and so 
𝙸_𝐏^-1𝙸𝙸_𝐏=𝙸𝙸_𝐏=(κ_1,κ_2) is exactly the Hessian matrix of F.
For some ϵ>0  small enough, consider on  the two planes  parallel to T_𝐏(𝒟) at distances  ±ϵ from T_𝐏(𝒟), 
the curves 𝒞_±={ (x̃,ỹ),  𝐅(x̃,ỹ)∈ T_𝐏(𝒟) | F(x̃,ỹ)=±ϵ}. 
It can be shown <cit.> that   the first-order approximation of the  intersections of 𝒟  with the two parallel planes is 
the union of two conics  (one real and one virtual) with equations 𝙸𝙸_𝐏(x̃,ỹ)=±2ϵ.
This union is known as the Dupin indicatrix when written in canonical form (i.e., by replacing 2ϵ by 1).






The real Dupin conic characterizes the local shape of 𝒟 and gives local information on the first-order geometry of the surface, 
at least at  points where the conic is non degenerate. It specializes to a parabola if the Gauss curvature vanishes i.e., κ_1κ_2=0, 
to an ellipse if κ_1κ_2>0, and  to an hyperbola if κ_1κ_2<0, see Fig. <ref>. 
Points are said to be elliptic, hyperbolic or parabolic respectively.




Let us now focus on the  Dupin central conics i.e., the real ellipse and real hyperbola.
We do not consider the case of a parabola since, in real experiments, the condition κ_2=0 will never be verified exactly.

Various measures can describe such a  conic shape <cit.>.
We introduce the CS quantity (<ref>)  that provides a  unified way of treating ellipses and  hyperbolas using the following nice interpretation.

Let the semi-major and semi-minor axes of the Dupin  central conic be r_maj=max(ρ_1 ^2, ρ_2 ^2) 
and r_min=min(ρ_1 ^2, ρ_2 ^2) respectively, where   ρ_1 and ρ_2 are the radii of curvature of the curves obtained 
through the  normal section of 𝒟 at 𝐏 along the principal directions.


CS in (<ref>)  is the squared ratio between the eccentricity e of the Dupin
 conic and  its minor axis (due to lack of space, the straightforward proof is omitted): 




    κ_1-κ_2=(e/r_min)^2

where

    e=√(
    1 ± r_min^2/r_maj^2) with
    -: if the conic is an ellipse
    
    +: if the conic is an hyperbola
 
 









The eccentricity can  be interpreted as the fraction of the distance along the semimajor axis at which the focus lies.







The CS quantity is normalized by additionally dividing the squared eccentricity  by the squared semi-minor axis.
Note that (<ref>) works for Dupin ellipses as well as Dupin hyperbolas.
The  curvilinear saliency CS  is large when κ_1 ≫κ_2, that is in presence of distant foci and so a highly elongated ellipse or a “squashed”
hyperbola, see Fig. <ref>.  
This occurs e.g., when the point is located on a depth “discontinuity”.
In turn, when κ_1 ≃κ_2, the conic approaches a circle and the distance between foci becomes very small.

A simple way to compute the curvilinear  feature:   
After  algebraic manipulations, it can be  shown   that

𝙸_𝐏^-1𝙸𝙸_𝐏=1/α𝖬
where


    𝖬≜[ (Z^2_y  + 1) Z_xx - Z_x Z_y Z_xy  (Z^2_y  + 1)Z_xy - Z_x Z_y Z_yy; (Z^2_x  + 1) Z_xy - Z_x Z_y Z_xx  (Z^2_x  + 1)Z_yy - Z_x Z_y Z_xy ]



 The squared curvilinear feature can be  computed as

    *CS^2   ≜‖∇_Z ‖^2(
    (𝖬)^2-4𝖬)
       =4‖∇_Z ‖^2(κ̅^2-K)

where κ̅ is the mean curvature of  𝒟 and K its
Gaussian curvature.


Let
λ_1 and λ_2 be the ordered eigenvalues (λ_1 ≥λ_2,)
computed from -𝖬 as defined in (<ref>),
to have κ_i=‖∇_Z ‖λ_i, i=1,2.
As the two eigenvalues of -𝖬 write

    λ_±= 1/2(-(𝖬) ±√( (𝖬)^2-4𝖬))

 we have λ_1=λ_+ and λ_2=λ_-. Since  as well
as

    κ_1-κ_2=‖∇_Z ‖(λ_1-λ_2)

the  squared curvilinear saliency is then defined.



However, the rely on the highest or smallest principal
curvature alone is not adequate
for defining accurate ridges <cit.>.
In Fig. <ref>, we show the different detections obtained
using the minimum or the maximum principal
curvature, as proposed by <cit.>.
The maximum principal
curvature provides a high response only for dark lines on a
light background, while the minimum gives the higher answers for the light
lines on a dark background.
The difference of the principal
curvatures, κ_1 - κ_2, improves robustness
as it responds in both settings.





§ IMAGE REPRESENTATION








 §.§ Proposed Curvilinear Features for Images


Let I(x,y) denote the value of the  image intensity function I:U⊂ℝ^2→ℝ at   image point (x,y).
Similar to the work of <cit.>, we assume the Lambertian shading model (<ref>). 
Let the  intensity image be treated as an intensity surface ℐ defined by the vector function

    𝐈(x,y)=[ x,y,I(x,y) ]^⊤

Remind that the unit normal is 𝐍(x,y)=𝐍̅(x,y)/‖𝐍̅(x,y)‖ where 𝐍̅ is defined in (<ref>),
and so, only depends on the depth Z(x,y) and its derivatives up to order-1. 









We now want to detect features in the intensity surface ℐ and check whether they are good candidates to be matched to curvilinear features detected 
in the depth surface  𝒟, w.r.t. a given camera pose. 
The key issue here is that detected features in ℐ can be matched to features detected in 𝒟 on the condition that 
both are based on measurements with the same order of derivation in  Z(x,y),  in order to yield a “compatible” matching that ensures  repeatablility. 
The fact that I  depends on   Z(x,y) and its derivatives up to order-1,  entails that the detection of  features in  ℐ must rely
on  order-1 variations of the surface 𝐈(x,y), e.g., on its differential along some adequate direction.

Consider  a point 𝐐=𝐈(x,y) on the image surface. Let   d𝐈_𝐐:U→ℝ^3 be the differential of 𝐈 at 𝐐. Given a unit direction   𝐯=[a,b]^⊤ in the image xy-plane, we have d𝐈_𝐐(𝐯)
=a𝐈_x+b𝐈_y=𝙹_𝐈𝐯

















where 𝙹_𝐈







is the Jacobian matrix of 𝐈 and  
𝐈_x=[1,0,I_x]^⊤ and 
𝐈_y=[0,1,I_y]^⊤, where

    I_⋆=1/2(𝐍·𝐬-1 )(𝐍_⋆·𝐬)

⋆ standing for either x and or y.
It is an order-1 measurement of the image surface variation  at 𝐐 and is compatible with  our curvilinear measurements of the depth surface 
(i.e., with same order of the derivatives of Z).

In order to a get a scalar measurement, we define the unit vectors
𝐓_1=𝙹_𝐈∇_I/‖∇_I‖
and 𝐓_2 by rotating 𝐓_1 by π/2. For α=1,2, we also define
    |μ_α|=‖d𝐈_𝐐(𝐓_α)
    ‖




















which is the differential of 𝐈  along unit direction 𝐓_α in the image plane.  
It can be easily seen that ∇_I/‖∇_I‖ is the eigenvector of 

    𝙹_𝐈^⊤𝙹_𝐈
    =
    [ 𝐈_x·𝐈_x 𝐈_x·𝐈_y; 𝐈_x·𝐈_y 𝐈_y·𝐈_y;         ]   =
    [ 1+(I_x)^2   I_x I_y;   I_x I_y 1+(I_y)^2;           ]
       =𝙸 +∇_I∇_I^⊤

associated with the largest eigenvalue μ_α.
It is worthy to note that the similarity between
the expression of the principal curvature computed for the depth surface, cf. (<ref>) and the formulae (<ref>).
Also, note that the matrix (<ref>) is that of  the first fundamental
form of  ℐ.
Clearly, the maximum and minimum values of the quadratic form ‖d𝐈_𝐐(𝐯)‖^2
correspond to the two  eigenvalues of the first fundamental form matrix given in (<ref>).

By a similar approach to <ref>, we can propose as feature the difference μ_1-μ_2
where μ_1≥μ_2.



























Let μ_1,μ_2 be the two eigenvalues of the first fundamental form matrix
𝙹_𝐈^⊤𝙹_𝐈 of ℐ,  ordered in descending
order.
Then, we have

    μ_1-μ_2=‖∇_ I‖^2



 We can deduce the ordered eigenvalues of 
𝙸^ℐ_𝐏
=𝙹_𝐈^⊤𝙹_𝐈
from those of   ∇_I∇_I^⊤, i.e., ‖∇_
I‖^2 and 0, so μ_1=‖∇_ I‖^2+1
and μ_2=1. Which ends the proof.


  Again, as in <ref>, we can describe the local shape of ℐ at 𝐐
  by means of the  eccentricity of a conic, here given by
 the quadratic form 𝐯^⊤𝙹_𝐈^⊤𝙹_𝐈𝐯=±
1. How can we interpret this conic? The first order Taylor expansion for
infinitesimal changes (dx,dy) in the vinicity of
  𝐐=𝐈(x,y)
 yields
 
    𝐈(x+ d x, y+dy) -𝐈(x,y) ≈𝙹_𝐈
    [d
    x,d y ]^⊤

For any unit direction 𝐯=[a,b]^⊤
in the xy-plane, the quadratic form 𝐯^⊤𝙹_𝐈^⊤𝙹_𝐈𝐯
returns the  linear part g of growth in arc length
from 𝐈(x,y) to 𝐈(x+a,y+b). Therefore, we have

    g^2=‖d𝐈_𝐐((dx,dy)‖^2
    =𝐯^⊤𝙹_𝐈^⊤𝙹_𝐈𝐯
 















































































































 




































































































































An important remark that we highlight here and not mentioned in <cit.> is the following one. 
The AVG feature defined in (<ref>) is actually the Frobenius norm of the Jacobian matrix 𝙹_𝐍 of the map 𝐍(x,y), see (<ref>),
up to  constant γ.  
Clearly, this describes the second-order behaviour of the surface  ℳ relative to the  normal at one of its points in the immediate vicinity of this point. 
Using the result in (<ref>), (<ref>) and (<ref>), we can claim that the extracted feature in the virtual image only depends on X, Y,Z 
and their derivatives up to order-2. 
It is consistent (regarding the considered orders of the derivatives of X,Y,Z) with the feature ‖∇_ I‖ = √(I_x^2+I_y^2) 
detected in the intensity image, where I_⋆, with ⋆∈{1,2} is given in (<ref>). 









































































































 §.§ Multi-Curvilinear Saliency (MCS)



Multi-scale helps to detect important structures as well as small details.
In consequence, in this paper, we compute the curvilinear saliency images in a multi-scale space.
To build the scale pyramid, an edge-preserving smoothing approach, named anisotropic diffusion filter <cit.>, is used in order to avoid oversmoothing.
In fact, this filter tries to separate the low frequency components (i.e, sharp edges) from the high frequency components (i.e., textures) 
by preserving the largest edges in an image.

Contrary to depth images which represent textureless 3D shapes, intensity images are composed of shape and texture components.
Consequently, the curvilinear saliency (CS) estimated from intensity images is affected by the textured regions.
Our idea is to put forward the assumption that multi-scale analysis can discriminate between keypoints 
(those with high CS value in the image) due to shape and keypoints due to texture.
At a coarse level, edges detected are reliable but with a poor localization and they miss small details.
At a fine level, details are preserved, but detection suffers greatly from clutters in textured regions.
In addition, the CS 
values of small details and textures are high in the coarse level,
whereas these values become lower in the finest levels.
To combine the strengths of each scale, the CS 
value of each pixel over n scales is analyzed.
If this value in all scales is higher than a threshold T, the maximum curvilinear saliency (*MCS) value of this pixel over all scales is then kept.
This threshold is a function of the number of the smoothed images, n, (i.e., T = e^-n: when n is small, then T is a big value and vice versa).
However, if the CS 
value is lower than T in one level, it is considered as a point that belongs to a texture (or a small detail) point,
thus it is removed from the final multi-scale curvilinear saliency, *MCS, image.


Adding this multiscale step should help to reduce the impact of the texture on the point of interest detection. 
However, in the next section, we propose to introduce the principle used for estimating focus map 
in order to increase the robustness to the background and to the presence of the texture. 



§ ROBUSTNESS TO BACKGROUND AND TO TEXTURE


Before introducing the proposed improvement, we briefly present existing works about texture detection and, in particular, about focus curve estimation. 


 §.§ Extraction of texture: state of the art


Various methods, such as <cit.> have been proposed for extracting the texture from a natural image.
In these approaches, a given image is separated into two components while  preserving edges. In <cit.>, Aujol et al. proposed a variational model based on total-variation (TV) energy for extracting the structural part. 
In <cit.>, the authors proposed an algorithm in the field of scale space theory. 
This technique is a rolling guidance method based on an associated depth image to automatically refine the scale of the filtering in order to preserve edges. 
A structure-preserving image smoothing approach is introduced in <cit.>.
This method locally analyzes the second order feature statistics of a patch around a pixel.
The algorithm used a 7-dimensional feature vector that contains intensity, orientation and pixel coordinates.

Finally, under a condition that the images contain smooth illumination changes and Lambertian surfaces,
<cit.> proposed an intrinsic image decomposition model which explicitly determines a separate texture layer, as well as the shading layer and the reflectance layer.
The method is based on surface normal vectors generated from an RGB-D image.
All these works first smooth the intensity image as a pre-processing stage and then extracting the shape from that image relying on  prior knowledge.
And, to sum up, most of these methods for structure-texture decomposition are analogous to the classical signal processing low pass-high pass filter decomposition.
However, even if it is correct to consider that the structure part of an image contains strong edges, the texture can also contain medium and high frequencies and the texture can only be partially removed.
Another possibility is to consider focusness, related to the degree of focus.

Usually focusness is defined as inversely proportional to the degree of blur (blurriness) <cit.>.
It is a very valuable tool for depth recovery <cit.> but also for blur magnification, or for image quality assessment.
Blur is usually measured in regions containing edges, since edges would appear in images as blurred luminance transitions of unknown blur scale <cit.>. Then the estimation of the blur can be propagated to the rest of image.
Since blur occurs for many different causes, this task is challenging and, in the literature, many methods on focus map estimation have been proposed.
In <cit.>, the authors identify blur as focal blur, induced by  finite depth of field, as penumbra blur or shading blur and then estimate the blur scale. 
In the context of matting method <cit.>, the blur ratio for every pixel corresponds to the ratio between the gradients of the input image and the re-blurred images.
In <cit.>, the authors use the K nearest neighbors (KNN) matting interpolation under the assumption that depth is locally uniform.
However, blurring can also appear with edges caused by shadows and glossy highlights can also produce error in focus estimation.


To remove errors induced by these other sources of blur, <cit.>  used a cross bilateral filtering and estimation of sharpness bias and <cit.> used multi-scale.
In the field of saliency detection, another way to estimate the amount of blur consists in computing the defocus blur between two Difference of Gaussian (DoG) images 
in multi-scale levels <cit.>.

An optimization algorithm is used by minimizing the difference between the blurriness of a pixel and the weighted average of blurriness of neighboring pixels.
In a Markov Random Field (MRF) formulation, a local contrast prior based on comparing local contrast and local gradient is also introduced in <cit.>.
In <cit.>, the authors propose to use the ratio between principle components, that is a weighted mixture of the spectral components.
Moreover, the weights are proportional to the energy in the spectral component.
Some algorithms also use the analysis of localized Fourier (Gabor filtering) spectrum <cit.>.
In addition, smoothness constraints and image color edge information are taken into account to generate a defocus map for trying to preserve discontinuities 
on the transitions between objects.



Following all these aforementioned approaches, we can find that most of the existing algorithms <cit.>
depend on measuring the blur amount using the ratio between the edges in two different scale levels (i.e., the original image and the re-blurred image).
In consequence, we propose to use the ratio between the two curvilinear saliency images that contain robust edges in different scales to determine the blur amount based on the methods developed in <cit.>. 
For the multi-scale aspect, our approach is inspired by the principles explained in <cit.>.




 §.§ Removing background with focus curves: state of the art


Based on the mapping between the depth of a point light source and the focus level of its image,
Shape From Defocus (SFD) approaches recover the 3D shape of a scene from focused images that represent the focus level of each point in the scene <cit.>.
We can also notice that the focus (defocus) maps can be also used as an alternative for depth map, like in existing Adobe tools <cit.>.
Consequently, it seems interesting to introduce what we call the detection of “focus curves” that capture blurriness in images.
More precisely, focus curves mean that we estimate the scale of blur at the curvilinear saliency feature of the original image and we suppose that these features should be only related to discontinuities.

Focal blur occurs when a point is out of focus, as illustrated in Fig <ref>.
When the point is at the focus distance, d_f, from the lens, all the rays from it converge to a sharp single sensor point.
Otherwise, when d ≠ d_f, these rays generate a blurred region in the sensor area.
The blur pattern generated by this way is called the circle of confusion (CoC) whose diameter is denoted c.






In <cit.>, the defocus blur can be modeled as a convolution of a sharp image with the point spread function (PSF) as shown in Fig. <ref>.
The PSF is usually approximated by a Gaussian function g(x, σ), where the standard deviation σ∝ c measures the blur amount and is proportional to the diameter of the CoC:

    c = |d-d_f|/df/d-f,

where d, d_f, f are focus distance, defocus distance and focal length respectively as shown in Fig. <ref>.
A blurred edge i(x) is then given by

    i(x)=f(x)⊗ g(x,σ)

where f(x)=Au(x)+B is an ideal edge where u(x) is the step function. 
The terms A and B corresponds to the amplitude and the offset of the edge, respectively. Note that the edge is located at x=0.

In <cit.>, the blur estimation method was described for 1D case. The gradient of the re-blurred edge is:

    ∇ i_1(x) =   ∇ (i(x)⊗ g(x,σ_0))
    
    =    ∇ ((Au(x)+B) ⊗ g(x,σ) ⊗ g(x,σ_0))
    
    =    A/√(2π (σ + σ_0))exp(-x^2/2(σ^2 + σ_0^2))

where σ_0 is the standard deviation of the re-blur Gaussian kernel. 
Thus, the gradient magnitude ratio between the original and the re-blurred edges is:

    R    = |∇ i(x)|/|∇ i_1(x)|
       = √(σ^2 + σ_0^2/σ^2)exp -(x^2/2(σ^2)-x^2/2(σ^2 + σ_0^2)).


It can be proved that the ratio is maximum at the edge location (x=0) and the maximum value is given by:


    R = √(σ^2 + σ_0^2/σ^2)


Finally, given the maximum value R at the edge locations, the unknown blur amount s can be calculated using:


    s = σ_0/√(R^2-1)




 §.§ Focus curves based on Curvilinear Saliency – Multi Focus Curves (MFC)


We suppose that using focus can help to remove the background and using multiscale can help to reduce the influence of the texture in the same way as in section <ref>.
So, we propose to use the curvilinear saliency computation instead of the edge response to estimate the focus curves of an input image.
In addition, we propose to estimate focus curves in multi-scales rather than in one scale as proposed in <cit.> to have scale invariant focus salient curves.
In addition, we combine all information gotten from different blurring scales.

Assume the original pixel in an image is blurred as F(x,y)= I(x,y)⊗ g(x,y,σ). Thus to get  the curvilinear saliency, the structure tensor is calculated as:

    𝖲_𝖳 =    f(∇((I(x,y) ⊗ g(x,y,σ))))
    
        =   ( [  I_x^2 I_xI_y; I_xI_y  I_y^2 ])⊗ g(x,y,σ)


If the Hessian matrix is expressed with eigenvectors and eigenvalues, we obtained:

    𝖧 =([ [ e_1 e_2 ]] [ [ λ_1   0;   0 λ_2 ]] [ [ e_1^T; e_2^T ]] ) ⊗ g(x,y,σ)


 
    𝖧 = ( (λ_1-λ_2) e_1e_1^T + λ_2(e_1e_1^T+e_2e_2^T))⊗ g(x,y,σ)

The curvilinear saliency can be described as:

    CS = (λ_1 -λ_2) ⊗ g(x,y,σ)

In particulary, the curvilinear saliency can be directly computed as:

    CS = α ((I_x^2+I_y^2)) ⊗ g(x,y,σ)

The re-blurred curvilinear saliency image, named CS_i, in multi-scale can then be defined as:

    CS_i = α ((I_x^2+I_y^2)) ⊗ g(x,y,σ) ⊗ g(x,y,σ_i), i=1,2,...,n

where n is the number of scales.

Consequently, the ratio between the original and re-blurred curvilinear saliency is:

    R_i = CS_i/CS = σ^2 + σ_i^2/σ^2exp -(x^2+y^2/2(σ^2)-x^2+y^2/2(σ^2 + σ_i^2))

Within the neighborhood of an pixel, the response reaches its maximum when x = 0 and y =0, thus:
 
    R_i|_0 = CS_i/CS = σ^2 + σ_i^2/σ^2 = 1+σ_i^2/σ^2

Finally, given the maximum value R_i in each scale level, the unknown blur amount s_i can be calculated using

    s_i = σ_i/√(R_i|_0-1),


For n scales, we compute n-1 focus curve scales by using the ratio between curvilinear saliency of the coarse level (i.e., the original image) and the next scale levels.



By following the same remarks as in section <ref>, we define 

Multi Focus Curves (MFC) that correspond to the fusion of all the focus curves into one map by keeping only the pixels that have
a high focus value in all the n-1 scales (i.e., a high value means a value bigger than T=e^-n, chosen in the same way as section <ref>).
If the pixel has a high value at all scales, the maximum value of the scale of blur is taken into account to build the final multi-scale curve map :

    MFC = 1/max_i (s_i).


In conclusion, the highest values of the estimated MFC indicate edges that have low blurring (i.e., sharp edges).
On the contrary, low values indicate ones that have high blurring.
Consequently, we expect that focus curves highlight salient curvilinear saliency in images that are approximately similar to the detected curvilinear saliency features in depth images.




§ EXPERIMENTS FOR FEATURE DETECTION



 §.§ Comparison with existing methods


One of our most important objectives in this work was to introduce a detector that is more repeatable between 2D images and 3D models than classical detectors in the litterature. 
In consequence, we compare the features detected on 3D models with the proposed curvilinear saliency detector with features detected on real images with these three 2D detectors:
Image Gradient (IG), Multi-scale Curvilinear Saliency (MCS) and  Multi-scale Focus Curves (MFC). 
In addition, we measure the repeatability between the two others 3D model detectors, i.e. Average Shading Gradient (ASG) <cit.> and Hessian Frobenius Norm (HFN), 
and the same three 2D detectors. And then, we compared MFC and MCS with nine classical 2D detectors:

 
  * Edge detectors:  (i) Sobel,
  (ii) Laplacian of Gaussian (Log),
  (iii) Canny <cit.> and
  (iv) Fuzzy logic technique <cit.>;
  
  * Corner detectors:  (v) Harris detector based on auto-correlation  analysis and
 (vi) Minimum Eigenvalues detector based on analysis of the Hessian matrix <cit.>;
 
  * Multi-scale detectors: (vii) SIFT, Scale Invariant Feature Transform <cit.>, that uses the analysis of difference of Gaussian,
 (viii) SURF, Speeded Up Robust Features, a multi-scale technique based on the Hessian matrix <cit.> and
  (ix) a multi-scale Principal Curvature Image (PCI) detector <cit.>.





 §.§ Evaluation criteria

The eleven 2D detectors are evaluated with two measures:

 
  * Intersection percentage (IP): the probability that a 2D intensity-based key feature can be found close to those extracted in a depth image <cit.>.
 
  * Hausdorff Distance (HD): the classical measurement is defined for two point sets A and B by: 
 
    HD(A,B)=max(h(A,B), h(B,A)),

 where h(A,B) = max_a∈ Amin_b∈ B∥ a-b ∥. The lowest the distance, the most similar the two sets.


 

 §.§ Datasets


Two datasets are evaluated:

 
  * Web collection: we have collected 10 objects and 15 real images of each object on the web by choosing views as close as possible to the views used for the generation of the depth images.
Moreover, to highlight the robustness of the approach to different acquisition conditions, many real images of a similar model are taken.

  * PASCAL3D+ dataset <cit.>: it is used in order to assess scalability. it contains real images corresponding to 12 rigid objects categories.
We have computed average results for all non occluded objects in each category, i.e. around 1000 real images and 3 or more reference models per category.
The real images are acquired under different acquisition conditions (e.g., lighting, complex background, low contrast).
We have rendered the depth images of the corresponding 3D CAD model using the viewpoint information from the dataset.
Only non-occluded and non-truncated objects in the real images were used. Furthermore, we choose 3D textureless objects (available online [<http://tf3dm.com/>]),

For all the tested 3D models, we have rendered depth images using MATLAB 3D Model Renderer [<http://www.openu.ac.il/home/hassner/projects/poses/>].



 §.§ Analysis of the results



As shown in tables <ref> and <ref>, and as expected, the proposed approach using focus curves based on curvilinear saliency, named MFC, 
is able to find the highest number of features in the intersection with the features detected on real images captured under different textures and lighting conditions.
More precisely, MFC obtains an average mean intersection percentage greater than 56% whereas for MCS and PCI, it is respectively greater than 50% and 44%, 
for the web collection dataset.
With the PASCAL+3D dataset, MFC also yields the highest mean average IP among all the tested detectors that is 46%.






In addition, as shown in <ref> and <ref>, the average Hausdorff Distance (HD) with MFC is less than 24 and with MCS, less than 32.

On the contrary, the other detectors do not reach high repeatability scores. 

All these quantitative results support that MFC is able to detect curvilinear saliency features that are more repeatable 
between an intensity image and its corresponding depth image than the state of the art.







In the rest of this section, we illustrate the results for the most significant dataset, the PASCAL3D+, in order to avoid redundancy in the explanations. 
In Fig. <ref>, two examples of results obtained for the PASCAL3D+ dataset <cit.> are given.
More precisely, the repeatability percentage between the three comparable 3D detectors, i.e. MFC, MSC and Image Gradient (IG), and 
the three comparable 2D detectors, Hessian Frobenius Norm, Average Shading Gradient and 
CS, is presented. 
These results highlight that image gradients are effected by texture.
Moreover, MCS improves the repeatability between depth images and real images, compared to IG.
And, as expected, MFC still yields the best repeatability scores. 
Among the detectors used for depth images, the Curvilinear Saliency detector yields the best repeatability scores between the three intensity-based 2D detectors. 
In conclusion, using CS with MFC gives the best repeatability among all the other possible combinations. 

In Fig. <ref>, we show some visual results with the PASCAL+3D dataset. 





The MFC and MCS detectors were applied on real images and their corresponding depth images. 
As shown, MFC can reduce a lot of edges belonging to texture information and can provide an approximation of the object shapes present on depth images.





§ REGISTRATION OF 2D IMAGES TO 3D MODELS


In this section, we register a 2D query image to a 3D model by finding the closest view d between all the rendered images of the 3D model 
d_k, k=1… N, where N the number of rendered views (i.e., depth images). 
We suppose that the object to recognize is contained in a bounding box and we want to estimate the 3D pose. 
Estimating the pose consists in estimating the elevation and the azimuth angles, respectively (h)  and (a), and the distance between the model and the camera, (v). 
For each 3D model, we generate depth images from near uniformly distributed viewing angles around a sphere by changing h, a and v to have N views per model. 
The choices for these terms are explained in paragraph <ref>.


To describe an object in a photograph and in all the rendered depth images, we naturally expand the famous classical HOG, Histogram Of Gradient, 
widely used in the literature <cit.>, to work on curvilinear saliency. 
HOG is  used in a sliding window fashion in an image to generate dense features based on binning the gradient orientation over a spatial region. 
Indeed, both in rendered depth images and in photographs, the orientation of the curvature and the magnitude of the curvilinear saliency is used for building the descriptors. 
For depth image representation, we multiply CS by the eigenvector e_H_1 corresponding to the largest eigenvalue of the matrix 𝖬 in (<ref>):

    CS = CS.e_H_1.

In turn, for photographs, for MCS, we multiply MCS values by the eigenvector e_S_1 corresponding to the curvilinear saliency λ_1-λ_2, 
as shown in (<ref>):

    MCS = MCS.e_H_1.

In addition, for MFC, we also multiply MFC values with the eigenvector e_S_1:

    MFC = MFC.e_S_1.

Using the same principle of HOG presented in <cit.>, 
we propose a descriptor that contains the orientation of the curvature and the value of curvilinear saliency and the magnitude of CS, MCS and MFC 
for an image are binned into sparse per-pixel histograms.

Given the HOG descriptor from a 2D query image D_q, we compute the HOG descriptors of the rendered images D_d_N, with N rendering depth images. 
In order to compare D_q to every D_d_N, the similarity scores are computed as proposed in <cit.>:

    𝐒_hog(k,h,a,v) = (𝐃_d_N - μ_s)^TΣ^-1𝐃_𝐪,

where, k=1...N, Σ and μ_s are, respectively, the covariance matrix and the mean over all descriptors of the rendered images.
For the registration process, evaluating 𝐒_hog(k,h,a,v) can be done by computing the probability of the inverse of the inner product between D_q and a transformed set of descriptors. The 𝐒_hog(k,h,a,v) probability is then maximized in order to find the closest corresponding views of the query image.

We also evaluate a global similarity by measuring how well each individual detected point in an image is able to be matched with a corresponding detected point in the depth map, i.e., how well each image detected points are repeatable.
More precisely, this repeatability scores, Rep, normalized between 0 and 1, is the probability that key features in the intensity image are found close to those extracted in the depth image Rep_d_N  q.
Since the closest view should have a high repeatability scores in comparison to other views, the dissimilarity based on repeatability scores is defined by R_d_i = 1-Rep_d_i q. If we denote 𝐑_d_i the repeatability scores of N rendered views of a model and a given image, the similarity 𝐒_rep is defined by:


    𝐒_rep(k,h,a,v) = exp (-(𝐑_d_i-μ_r)^2/2  σ_r^2).

where μ_r is the mean value of 𝐑_d_N and σ_r is the standard deviation (i.e., in this work σ_r = 0.1).

Finally, by combining all HOG feature similarities and the similarity based on the repeatability, the probability of the final similarity is given by:

    𝐒(m,h,a,v) =  𝐒_hog(k,h,a,v)⊙𝐒_rep(k,h,a,v).

where ⊙ is the Hadamard product. Based on calculating 𝐒(k,h,a,v), we select at least the highest three correspondences to estimate the full pose. From the selected three views, the logically ordered or connected views (i.e., coherent views) are firstly selected.
We then find minimum and maximum values of h, a and v of the corresponding views.
Additional views are then generated in the vicinity of the selected views that is between the minimum and maximum values of the three parameters with small steps (e.g., δ h = 5^o, δ a = 5^o and δ v = 5 cm).

The process is again repeated for these ranges to find the closet view to the object in a query image until convergence. 

Assume the ground-truth transformation matrix (𝐓_𝐠) containing rotation (𝐑_𝐠) and translation (𝐭_𝐠) matrices:
𝐓_𝐠 = [[ 𝐑_𝐠 𝐭_𝐠;   0   1;     ]].
In addition, the estimated transformation matrix (𝐓_𝐞) containing rotation (𝐑_𝐞) and translation (𝐭_𝐞) matrices:
𝐓_𝐞 = [[ 𝐑_𝐞 𝐭_𝐞;   0   1;     ]].
Then the matrix 𝐌 is computed: 𝐌 = 𝐓_𝐞^-1×𝐓_𝐠. The error between the two transformation matrix is: E =  M_i - M_j, where i≠ j. Thus, the convergence criterion is E ≤ε, ε is a very small value (in this work, ε=0.05).



§ EXPERIMENTS FOR POSE ESTIMATION





 §.§ 3D models representation and alignment


Matching photographs and rendered depth images requires a complete 3D model representation.
Each depth image represents a 3D model from different viewpoints.
Actually, we need to have a large number of depth images to completely represent a 3D model.
However, this yields a massive execution time.
Consequently, we have orthographically rendered N depth images (around 700 in our experiments) from approximately uniformly distributed viewing angles h and a and the distance v
(i.e., in these experiments, h is increased by a step of 50^o, the azimuth angle, 20^o, and the distance, 0.3 m, for a range between 0 and 2 m).

In addition, we need a parametrization of the alignment of the model view in a depth image with the object detected in a color image.
Consequently, when we compare two models, we need to compute the optimal measure of similarity, over all possible poses.
To do this, each model is placed into a canonical coordinate frame, normalizing for translation and rotation.
Since we know the centroid of them, the models are normalized for translation by shifting them so that the center of mass is aligned with the origin.
Next, the two models are normalized for rotation by aligning the principal axes of the model with the x-, and y-axes.
It defines the ellipsoid that best fits the model. By rotating the two sets of points so that the major axis of the ellipsoid is aligned with the x-axis,
and the second major axis aligns with the y-axes, we obtain a model in a normalized coordinate frame.
We use Principle Components Analysis, PCA, to find the orientation of the major axis of the ellipse.
The set of points of the model is rotated by the difference of the direction of the two major axes.
After normalization, the two models are (near) optimally aligned and can be directly compared in their normalized poses.

In addition, in this paper, the HOG descriptor is quantized into 9 bins, exactly as proposed in <cit.>.
The photograph and each depth image is divided into a grid of square cells
(i.e., in this work, the image is divided into 8× 8[Different grids were tested: 4× 4, 8× 8 and 16× 16. The grid with 8× 8 size yields the best precision rate.]).
For each cell, the curvilinear saliency, focus curves or image gradients, histograms are aggregated by weighting them with their respective magnitude.



 §.§ Analysis of the results



 

For pose estimation or even for object recognition, it is needed that the probability that key features in the photograph are found close to those extracted in the depth image must be high when the photograph and the depth image come from the same viewpoint.

In the first experiment that is presented, we try to illustrate that the correct view provides the highest repeatability score, in particular with the proposed detector.
In consequence, the repeatability scores of features extracted on photograph and depth images with MCS, SIFT and Sobel edge detector, 
for the correct view and the other views around it, were calculated and
shown in Fig. <ref> with one example of the dataset.
As shown, the three detectors yield the highest repeatability score with the correct viewpoint (even if the difference between views is small, like with Sobel).
In addition, as expected, the score is gradually diminished whenever it is far from the correct viewpoint.
The most important result is that MCS gives the highest differences between the correct view and all the other views, so, 
it illustrates that it is the most adapted detector for pose estimation based on 2D/3D registration.
This result is quite coherent because SIFT was designed to be robust against many changes, many difficulties and
it induces that the differences should be lower than MCS that is designed to be efficient in the case of 2D/3D matching.


In addition, the second experiment was performed with the Pascal+3D dataset.
For each category of objects, we compute the precision rate for detecting the correct view after using 

the three aforementioned methods for 3D model representations, i.e., Curvilinear Saliency (CS), Average Shading Gradient (ASG) and Apparent Ridges (AR) <cit.>, 
against the three techniques for intensity image representation, i.e., Image Gradient(IG), Multi Curvilinear Saliency (MCS) and Multi Focus Curves (MFC).
As shown in table <ref>, the registration between our Curvilinear Saliency (CS) representation of the 3D model and the multiscale focus curves
(MFC) extracted on corresponding images outperforms all other variations of the tested methods.
This confirms that curvilinear saliency representation computed from the depth images of a 3D model can capture the discontinuities of surfaces.
In addition, MFC can reduce the influence of texture and background components by extracting the edges related to the object shape in intensity images rather than MCS.

Furthermore, the precision rate is reduced by more than 25% compare with ASG and IG to represent 3D models and intensity any other image representation.
Apparent ridge rendering yields the smallest registration accuracy with the three image representations among all the 3D model representation techniques.
However, using ASG with untextured 3D models against MFC and MCS also increases correct pose estimation rate.
All these results  indicate that average shading gradients computed from the normal map of an untextured geometry is a good rendering technique for the untextured geometry.
However, image gradients are not the good representation of intensity images in order to match with rendering images because it is affected by image textures. 
All these results are confirmed by the third experiment, in table <ref>, where the details are given about the precision of the pose estimation 
in terms of elevation, azimuth, yaw angles and distance. 





In the next experiment, in Fig. <ref>, we show the precision of the registration of images among the top r similarities, i.e. 
we sort all the similarity scores obtained for all the views, and we analyse the r first highest similarities (more precisely, the 1, 3, 5, 10 and 20 first ranks). 

The correct pose is searched within this set of views.
As shown, the precision rate increased when the number of views is increased, for any combination of 3D model reprensentation and image representation. 
However, MFC yields the highest precision with the three tested methods of representing 3D models (i.e., CS, ASG, and AR).
In addition, MCS yields good precision values.
In fact, IG yields the smallest precision values because the edges detected with texture information have a bad influence on estimating the successful registration.



Finally, the Fig. <ref> shows some examples of the Pascal+3D dataset of correct registrations with the top-ranked pose estimation.
It can be seen that our system is able to register an image with a great variety of texture and viewing angle.
In addition, the proposed algorithm can register images regardless of light changes in images. 







§ CONCLUSION AND PERSPECTIVE


After an analysis of existing tools for 2D/3D registration, the major goal of this paper was to propose an approach for 2D/3D matching more adapted, and in particular more justified, 
than existing approaches. For that purpose, we also proposed an evaluation protocol based on repeatibility study. 
More precisely, for doing this matching, we have studied these two important aspects: how to represent the studied data, in 2D and in 3D, and then, how to compare them. 
In this context, we introduce a 3D detector based on curvilinear saliency and a 2D detector based on the same principle but adapted in multi-scale and combined with the principle 
of focus curves. 
The interest of this new method was also illustrated by quantitative evaluation on pose estimation and 2D/3D registration. All the results are very encouraging and the next 
step of this work is to use this registration to identify defaults on objects. For that purpose, we need to study the robustness of this work to missing part of objects and 
to adapt the registration process in consequence. 

1
Markelj2012 P. Markelj and D. Tomaževič and B. Likar and F. Pernuš. A review of 3D/2D registration methods for image-guided interventions. Medical Image Analysis, 16(3):642–661, 2012. 
Pomerleau2015 F. Pomerleau and F. Colas and S., Roland. A Review of Point Cloud Registration Algorithms for Mobile Robotics. Found. Trends Robot, 4(1):1–104, 2015. 
Krizhevsky2012 A. Krizhevsky and I. Sutskever and G. E. Hinton. ImageNet Classification with Deep Convolutional Neural Networks. In International Conference on Neural Information Processing Systems. 2012. 
Hao2015 H. Su and C. R. Qi and Y. Li and L. J. Guibas. Render for CNN: Viewpoint Estimation in Images Using CNNs Trained With Rendered 3D Model Views. In IEEE International Conference on Computer Vision. 2015. 
Dellepiane2008 M. Dellepiane and M. Callieri and F. Ponchio and R. Scopigno. Mapping highly detailed color information on extremely dense 3D models: the case of David's restoration, In Eurographics Conference – Cultural Heritage Papers. 
Wu2008 C. Wu and B. Clipp and X. Li and J. M. Frahm and M. Pollefeys. 3D model matching with Viewpoint-Invariant Patches (VIP). In IEEE Conference on Computer Vision and Pattern Recognition. 2008. 
Agarwal2009 S. Agarwal and Y. Furukawa and N. Snavely and I. Simon and B. Curless and S. Seitz and R. Szeliski. Building Rome in a Day. Communications of the ACM, 54(10):105–112, 2011. 
Campbell2001 R. J. Campbell and P. J. Flynn. A Survey Of Free-Form Object Representation and Recognition Techniques. Computer Vision and Image Understanding, 81(2):166–210, 2001. 
Choy2015 C. B. Choy and M. Stark and S. Corbett-Davies and S. Savarese. Enriching object detection with 2D-3D registration and continuous viewpoint estimation. In IEEE Conference on Computer Vision and Pattern Recognition. 2015. 
Plotz2015 T. Plötz and S. Roth. Automatic Registration of Images to Untextured Geometry Using Average Shading Gradients. International Journal of Computer Vision, 125(1):65–81, 2017. 
Irschara2009 A. Irschara and C. Zach and J. M. Frahm and H. Bischof. From structure-from-motion point clouds to fast location recognition. In IEEE International Conference on Computer Vision. 2009. 
Szeliski2002 D. Scharstein and R. Szeliski, Richard. A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms. International Journal of Computer Vision, 47(1):7–42,2002. 
Szeliski2010 R. Szeliski. Computer Vision – Algorithms and Applications. Springer, 2011. 
Dalal2005 N. Dalal and B. Triggs. Histograms of oriented gradients for human detection. In IEEE Conference on Computer Vision and Pattern Recognition. 2005. 
canny1986 J. Canny. A Computational Approach to Edge Detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8(6):679–698, 1986. 
harris88 C. Harris and M. Stephens. A combined corner and edge detector. In Alvey Vision Conference, 1988. 
Fischer2014 P. Fischer and T. Brox. Image Descriptors Based on Curvature Histograms. In German Conference on Pattern Recognition. 2014
Smith1997 S. Smith and J. Brady.  SUSAN - A New Approach to Low Level Image Processing. International Journal of Computer Vision, 23(1):45–78, 1997. 
Lowe2004 D. G. Lowe. Distinctive image features from scale-invariant keypoints. International Journal of Computer Vision, 60(2):91–110, 2004. 
Matas2002 J. Matas and O. Chum and U. Martin and T. Pajdla. Robust wide baseline stereo from maximally stable extremal regions. In British Machine Vision Conference. 2002. 
Tuytelaars2004 T. Tuytelaars and L. Van Gool. Matching Widely Separated Views Based on Affine Invariant Regions. International Journal of Computer Vision, 59(1):61–85, 2004. 
Mikolajczyk2004 K. Mikolajczyk and C. Schmid. Scale & Affine Invariant Interest Point Detectors. International Journal of Computer Vision, 60(1):63–86, 2004. 
Bay2008 H. Bay and A. Ess and T. Tuytelaars and L. Van Gool. Speeded-Up Robust Features (SURF). IEEE Computer Vision and Image Understanding, 110(3):346–359, 2008. 
Deng2007 H. Deng and W. Zhang and E. Mortensen and T. Dietterich and L. Shapiro. Principal Curvature-Based Region Detector for Object Recognition. IEEE Computer Vision and Image Understanding, 2007. 
Chang2005 C. H. Lee and A. Varshney and D. W. Jacobs. Mesh Saliency. ACM transactions on graphics, 24(3):659–666, 2005. 
Sattler2011 T. Sattler and B. Leibe and L. Kobbelt. Fast image-based localization using direct 2D-to-3D matching. In IEEE International Conference on Computer Vision. 2011. 
Ohtake2004 Y. Ohtake and A. Belyaev and H.-P. Seidel. Ridge-valley Lines on Meshes via Implicit Surface Fitting. ACM transactions on graphics, 23(3):609–612, 2004. 
Lee2007 Y. Lee and L. Markosian and S. Lee and J. F. Hughes. Line Drawings via Abstracted Shading. ACM transactions on graphics, 26(3), 2007. 
Park2013 Y. Y. Lee and K. Park and J. D. Yoo and K. H. Lee. Multi-scale Feature Matching Between 2D Image and 3D Model. In SIGGRAPH Asia. 2013. 
Judd2007 T. Judd, F. Durand, and E. H. Adelson. Apparent ridges for line drawing. ACM transactions on graphics, 26(3), 2007. 
Godil2011 A. Godil and A. I.  Wagan. Salient local 3D features for 3D shape retrieval. In SPIE. 2011. 
Paudel2014 D. P. Paudel and C. Demonceaux and A. Habed and P. Vasseur. Localization of 2D Cameras in a Known Environment Using Direct 2D-3D Registration. In IEEE Conference on Pattern Recognition, 2014. 
Xu2016 C. Xu and L. Zhang and L. Cheng and R. Koch. Pose Estimation from Line Correspondences: A Complete Analysis and a Series of Solutions. IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(6):1209–1222, 2017. 
Tamaazousti2011 M. Tamaazousti and V. Gay-Bellile and S. N. Collette and S. Bourgeois and M. Dhome. NonLinear refinement of structure from motion reconstruction by taking advantage of a partial knowledge of the environment. In IEEE Conference on Computer Vision and Pattern Recognition. 2011. 
Liu2005 L. Liu and I. Stamos. Automatic 3D to 2D registration for the photorealistic rendering of urban scenes. In  IEEE Conference on Computer Vision and Pattern Recognition. 2005. 
Ramalingam2009 S. Ramalingam and S. Bouaziz and P. Sturm and M. Brand. Geolocalization using skylines from omni-images. In IEEE International Conference on Computer Vision Workshops, 2009. 
Parlett1998 B. Parlett. The Symmetric Eigenvalue ProblemSociety for Industrial and Applied Mathematics (SIAM). 1998.
Gallier2001 J.-H. Gallier. Geometric Methods and Applications: For Computer Science and Engineering, (Texts in applied mathematics). Springer. 2001.
Clarkson2001 M. J. Clarkson and D. Rueckert and D. L. G. Hill and D. J. Hawkes. Using photo-consistency to register 2D optical images of the human face to a 3D surface model. IEEE Transactions on Pattern Analysis and Machine Intelligence, 23(11):1266–1280, 2001. 
Kiranpreet2010 K. Kiranpreet and V. Mutenja and I. S. Gill. Fuzzy Logic Based Image Edge Detection Algorithm in MATLAB. International Journal of Computer Applications, 1(22):55–58. 2010.
Lim2014 J. J. Lim and A. Khosla and A. Torralba. FPM: Fine Pose Parts-Based Model with 3D CAD Models. In European Conference on Computer Vision. 2014.
ours H. A. Rashwan and S. Chambon and P. Gurdjos and G. Morin and V. Charvillat. Towards Multi-Scale Feature Detection Repeatable over Intensity and Depth Images. In IEEE International Conference on Image Processing, 2016.
Elder1998 J. H. Elder and S. W. Zucker. Local scale control for edge detection and blur estimation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20(7):699–716, 1998.
Aubry2014 M. Aubry, D. Maturana, A. Efros, B. Russell and J. Sivic. Seeing 3D Chairs: Exemplar Part-Based 2D-3D Alignment Using A Large Dataset of CAD models. In IEEE Conference on Computer Vision and Pattern Recognition. 2014.
Yong2013 Y. Y. Lee, M. K. Park, J. D. Yoo, and K. H. Lee. Multi-Scale Feature Matching Between 2D Image and 3D Model. SIGGRAPH Asia, 2013.
Tomasi1994 J. Shi and C. Tomasi, Good features to track, in IEEE conference on Computer Vision and Pattern Recognition, 1994.
Paris2009 S. Paris and F. Durand. 2009. A Fast Approximation of the Bilateral Filter Using a Signal Processing Approach. International Journal of Computer Vision, 81(24–52), 2009.
Pentland1987 A. P. Pentland. A new sense for depth of field. IEEE Transactions on Pattern Analysis and Machine Intelligence 9(4):523–531, 1987.
Soonmin2007 B. Soonmin and F. Durand. Defocus magnification. Computer Graphics Forum, 26(3):571–579, 2007.
Zhuo2011 S. Zhuo and T. Sim. Defocus map estimation from a single image, Pattern Recongition, 44(9):1852–1858, 2011.
Tai2009 Y.-W. Tai and M. S. Brown. Single Image Defocous Map Estimation Using Local Contrast Prior. In IEEE International Conference on Image Processing. 2009.
Jiang2013 P. Jiang and H. Ling and J. Yu and J. Peng. Salient Region Detection by UFO: Uniqueness, Focusness and Objectness. In IEEE International Conference on Computer Vision, 2013.
Kumar2015 H. Kumar and S. Gupta and K. S. Venkatesh. Defocus Map Estimation from a Single Image using Principal Components. In International Conference on Signal Processing, Computing and Control. 2015.
Zhu2013 X. Zhu and S. Cohen and S. Schiller and P. Milanfar. Estimating Spatially Varying Defocus Blur From A Single Image. IEEE Transcations on Image Processing, 22(12):4879–4891, 2013. 
Wang2016 X. Zhang and R. Wang and X. Jiang and W. Wang and W. Gao. Spatially variant defocus blur map estimation and deblurring from a single image. Journal of Visual Communication and Image Representation 35:257–264, 2016.
Fischler1981 M. A. Fischler and R. C. Bolles. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography. Communications of the ACM, 24(6):381–395, 1981.
Aujol2006 J.-F. Aujol and G. Gilboa and T. Chan and S. Osher. Structure-Texture Image Decomposition—Modeling, Algorithms, and Parameter Selection. International Journal of Computer Vision, 67(1):111–136, 2006. 
Zhang2014 Q. Zhang, X. Shen, L. Xu and J. Jia. Rolling Guidance Filter. In European Conference on Computer Vision, 2014. 
Karacan2013 L. Karacan and E. Erdem and A. Erdem. Structure-Preserving Image Smoothing via Region Covariances. SIGGRAPH Asia, 2013. 
Jeon2014 J. Jeon and S. Cho and X. Tong and S. Lee. Intrinsic Image Decomposition Using Structure-Texture Separation and Surface Normals.  In European Conference on Computer Vision, 2014. 
PASCAL3D Y. Xiang, R. Mottaghi and S. Savarese. Beyond PASCAL: A Benchmark for 3D Object Detection in the Wild. In IEEE Winter Conference on Applications of Computer Vision (WACV), 2014.
doCarmo2016 M. P. do Carmo. Differential geometry of curves and surfaces. Prentice Hall, 1976. 
Hartley2004 R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision, Cambridge University Press, 2004. 
Koenderink1992 J. J. Koenderink and A. J. Van Doorn. Surface shape and curvature scales. Image and Vision Computing 10(8):557–564, 1992. 
Besl1992 P. J. Besl and N. D. McKay. A Method for Registration of 3-D Shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence. 14(2):239–256, 1992. 
Gray2006 A. Gray and E. Abbena and S. Salamon. Modern Differential Geometry of Curves and Surfaces with Mathematica. Chapman & Hall/CRC, 2006. 



——————————

By computing 𝐈_x=(1,0,I_x)^⊤ and 𝐈_y=(0,1,I_y)^⊤,
 it can be easily shown that
 
    𝙹_𝐈^⊤𝙹_𝐈
     =
     [ 𝐈_x·𝐈_x 𝐈_x·𝐈_y; 𝐈_x·𝐈_y 𝐈_y·𝐈_y;         ]   =
     [ 1+(I_x)^2   I_x I_y;   I_x I_y 1+(I_y)^2;           ]
       =𝙸 +∇_I∇_I^⊤

———————————

