2009 JPSS,Introduction of shape skewness parameter(s) in a probability distribution

BRIEF COMMUNICATIONRheological behaviour of ethylene glycol-titanate nanotube nanofluidsHaisheng Chen ÆYulong Ding ÆAlexei Lapkin ÆXiaolei FanReceived:11July 2008/Accepted:4February 2009/Published online:26February 2009ÓSpringer Science+Business Media B.V.2009Abstract Experimental work has been performed on the rheological behaviour of ethylene glycol based nanofluids containing titanate nanotubes over 20–60°C and a particle mass concentration of 0–8%.It is found that the nanofluids show shear-thinning behaviour particularly at particle concentrations in excess of *2%.Temperature imposes a very strong effect on the rheological behaviour of the nanofluids with higher temperatures giving stronger shear thinning.For a given particle concentration,there exists a certain shear rate below which the viscosity increases with increasing temperature,whereas the reverse occurs above such a shear rate.The normalised high-shear viscosity with respect to the base liquid viscosity,however,is independent of temperature.Further analyses suggest that the temperature effects are due to the shear-dependence of the relative contributions to the viscosity of the Brownian diffusion and convection.The analyses also suggest that a combination of particle aggregation and particle shape effects is the mechanism for the observed high-shear rheological behaviour,which is also supported by the thermal conductivity measure-ments and analyses.Keywords Rheological behaviour ÁEthylene glycol ÁTitanate nanotube ÁNanofluid ÁThermal conductivityNanofluids are dilute suspensions of particles with at least one dimension smaller than about 100nm (Choi 1995).Such a type of materials can be regarded as functionalized colloids with special requirements of a low-particle loading,a high-thermal performance,favourable flow/rheolgocial behaviour,and a great physical and chemical stability over a wide range of process and solution chemistry conditions.Nano-fluids have been shown to be able to enhance heat transfer (Choi 1995;Wang and Mujumdar.2007),mass transfer (Krishnamurthy et al.2006),and wetting and spreading (Wasan and Nikolov 2003),and have been a hot topic of research over the past decade (Wang and Mujumdar 2007;Keblinski et al.2005).Most published studies have focused on the heat transfer behaviour including thermal conduction (Choi 1995;Wang et al.1999;Wang and Mujumdar 2007;Keblinski et al.2005;Eastman et al.2001;He et al.2007;Ding et al.2006),phase change (boiling)heat transfer (Das et al.2003;Pak and Cho 1998),and convective heat transfer (Wang and Mujumdar 2007;Keblinski et al.2005;He et al.2007;Ding et al.2006,Chen et al.2008;Prasher et al.2006a and Yang et al.2005).Only few studies have been devoted to the rheological behaviour ofH.Chen ÁY.Ding (&)Institute of Particle Science and Engineering,University of Leeds,Leeds,UK e-mail:y.ding@pkin ÁX.FanDepartment of Chemical Engineering,University of Bath,Bath,UKJ Nanopart Res (2009)11:1513–1520DOI 10.1007/s11051-009-9599-9nanofluids(He et al.2007;Chen et al.2008;Prasher et al.2006a,b;Kwak and Kim2005;Lee et al.2006), although there is a large body of literature on suspensions rheology;see for example,Russel et al. (1991);Chow(1993);Petrie(1999),Larson(1999); Goodwin et al.(2000)l;Mohraz et al.(2004);Larson (2005);Egres and Wagner(2005);Abdulagatov and Azizov(2006).Particularly,there is little in the literature on the effect of temperature on the rheo-logical behaviour of nanofluids.Clearly,there is a gap in the current rheological literature for this type offluids.Furthermore,recent work has shown that the thermal behaviour of nanofluids correlates well with their rheological behaviour(Prasher et al.2006a, b;Chen et al.2007a;Abdulagatov and Azizov2006). In a recent study,we investigated systemically the rheological behaviour of ethylene glycol(EG)based spherical TiO2nanofluids(Chen et al.2007b).The results show that the nanofluids are Newtonian over a shear rate range of0.5–104s-1and the shear viscosity is a strong function of temperature,particle concentration and aggregation microstructure.This work is concerned about the rheological behaviour of EG based nanofluids containing titanate nanotubes (TNT).The specific objectives of the work are to investigate the effects of particle shape,particle concentration and temperature on nanofluids viscosity, and to understand the relationship between the rheo-logical behaviour and the effective thermal conductivity of nanofluids.It is for thefirst time that the rheological behaviour of a highly viscous EG based TNT nanofluids is investigated in a systematic manner.As will be seen later,the results of this work provide further evidence that the rheological measure-ments could provide information of particle structuring for predicting the effective thermal conductivity of nanofluids.The EG-TNT nanofluids used in this work were formulated by using the so-called two-step method with EG purchased from Alfa Aesar and TNT synthesized in our labs using a method described elsewhere(Bavykin et al.2004).The details of nanofluids formulation can be found elsewhere(Wen and Ding2005;He et al.2007;Chen et al.2007b). The TNT particles have a diameter(b)of*10nm and a length(L)of*100nm,giving an aspect ratio of(r=L/b)of*10.To avoid complications in interpreting the experimental results,no dispersants/ surfactants were used in the formulation.The nanofluids formulated were found stable for over 2months.The rheological behaviour of the nano-fluids was measured by using a Bolin CVO rheometer (Malvern Instruments,UK)over a shear rate range of 0.03–3,000s-1,a nanoparticle mass concentration of w=0–8%,and a temperature range of20–60°C (293–333K).The nanofluids were characterised for their size by using a Malvern Nanosizer(Malvern Instruments,UK)and a scanning electron microscope (SEM).The average effective particle diameter was found to be*260nm for all nanofluids formulated. This size is much larger than the equivalent diameter of the primary nanoparticles due to aggregation;see later for more discussion.Note that the particle size characterisation was performed both before and after the rheological measurements and no detectable changes to particle size were found.Figure1shows the viscosity of pure EG and EG-TNT nanofluids as a function of shear rate at 40°C.The results at other temperatures are similar.It can be seen that the EG-TNT nanofluids exhibit highly shear-thinning behaviour particularly when the TNT concentration exceeds*2%.Such behaviour is different from the observed Newtonian behaviour of EG-TiO2nanofluids containing spherical nanoparti-cles over similar shear rate range(Chen et al.2007b) where the base liquid,EG,is the same as that used in the current wok.The behaviour is similar to the observations of carbon nanotube nanofluids(Ding et al.2006)and CuO nanorod nanofluids(Kwak and Kim2005),although there are important differencesbetween them such as temperature dependence as will be discussed later.The shear-thinning behaviour of well-dispersed suspensions can be interpreted by the structuring of interacting particles(Doi and Edwards1978a,b and Larson1999).In a quiescent state,a rod-like particle has three types of motion due to Brownian diffusion: rotational(end-over-end)motion around the mid-point and translational motion in parallel or perpendicular to the long axis.For dilute suspensions with a number density,c,ranging between0and1/L3or volume fraction,u,ranging between0and1/r2),the average spacing between the particles is larger than the longest dimension of the rod,and zero shear viscosity can be approximated by gð0Þ%g0ð1þAÁcL3Þwith g0the base liquid viscosity and A,a numerical constant(Doi and Edwards1978a).For suspensions with 1/L3\c\1/bL2or1/r2\f/\1/r,the rod-like particles start to interact.The rotational motion is severely restricted,as well as the translational motion perpendicular to the long axis,and the zero shear viscosity can be estimated by gð0Þ%g0ð1þðBcL3Þ3Þ; with B a numerical constant(Doi and Edwards1978b). As a consequence,the zero shear viscosity can be much greater than the base liquid viscosity.The large viscosity is due to the rod-like shape effect and the viscosity is very sensitive to shear,which tends to align particles and hence the shear-thinning behaviour as shown in Fig.1.Note that the above mechanism can give a qualitative explanation for the experimental observations at low-shear rates and the shear-thinning behaviour as shown in Fig.1,it does not explain the high-shear viscosity of the nanofluids,which will be discussed later.It should also be noted that the criteria for classifying nanofluids given above need to be modified due to the presence of aggregates;see later for more discussion.Figure2shows the shear viscosity of4.0%EG-TNT nanofluids as a function of shear rate at different temperatures.The results under other concentrations are similar.It can be seen that the temperature has a very strong effect on the rheological behaviour of nanofluids with higher temperatures giving stronger shear thinning.For shear rates below*10s-1,the shear viscosity increases with increasing temperature, whereas the trend is reversed when the shear rate is above*10s-1.As mentioned above,this behaviour was not observed for carbon nanotube(Ding et al. 2006)and CuO nanorod(Kwak and Kim2005)nanofluids and we have not seen reports on such behaviour for nanofluids in the literature;see later for more discussion on the underlying mechanisms. Figure2also shows that the strongest shear thinning occurs at40–60°C,whereas very weak-shear thinning takes places at20–30°C.It is also noted that the shear viscosity of nanofluids at all temperatures investigated approaches a constant at high-shear rates.If the high-shear viscosity is plotted against temperature,Fig.3is obtained where the shear rate corresponding to the high-shear viscosity is taken as *2,000s-1.An inspection of all the data indicates that theyfit the following equation very well:ln g¼AþBÂ1000=TþCðÞð1Þwhere g is the shear viscosity(mPaÁs),T is the absolute temperature(K),and A,B and C areconstants given in Table1.Equation(1)takes a similar format as that widely used for liquid viscosity (Bird et al.2002)and for EG based nanofluids containing spherical particles(Chen et al.2007b).If the measured high-shear viscosity is normalized with respect to the shear viscosity of the base liquid, the relative increaseðg i¼ðgÀg0Þ=g0Þof the high-shear viscosity is found to be only a function of concentration but independent of temperature over the temperature range investigated in this work.The relative increments in the shear viscosities of nano-fluids containing0.5%,1.0%,2.0%,4.0%and8.0% particles are 3.30%,7.00%,16.22%,26.34%and 70.96%,respectively.Similar temperature indepen-dence of the shear viscosity was also observed for EG-TiO2and water-TiO2nanofluids containing spherical nanoparticles(Chen et al.2007b).The experimentally observed temperature depen-dence can be interpreted as follows.Given the base liquid and nanoparticles,the functional dependence of viscosity on shear rate is determined by the relative importance of the Brownian diffusion and convection effects.At temperatures below*30°C,the contribu-tion from the Brownian diffusion is weak due to high-base liquid viscosity.As a consequence,the shear dependence of the suspension is weak(Fig.2).The contribution from the Brownian diffusion becomes increasingly important with increasing temperature particularly above40°C due to the exponential dependence of the base liquid viscosity on temperature (Fig.3).At very high-shear rates,the Brownian diffusion plays a negligible role in comparison with the convective contribution and hence independent of the high-shear viscosity on the temperature.We now start to examine if the classical theories for the high-shear viscosity predict the experimental measurements(note that there is a lack of adequate theories for predicting the low shear viscosity).Figure4shows the shear viscosity increment as a function of nanoparticle volume concentration together with the predictions by the following Brenner &Condiff Equation for dilute suspensions containing large aspect ratio rod-like particles(Brenner and Condiff1974):g¼g01þg½ uþO u2ÀÁÀÁð2Þwhere the intrinsic viscosity,½g ;for high-shear rates has the following form(Goodwin and Hughes2000):½g ¼0:312rln2rÀ1:5þ2À0:5ln2rÀ1:5À1:872rð3ÞAlso included in Fig.4are the data for EG-TiO2 nanofluids with spherical nanoparticles(Chen et al. 2007b)and predictions by the Einstein Equation (Einstein1906,1911)for dilute non-interacting suspensions of spherical particles,g¼g01þ2:5uðÞ: It can be seen that both the Einstein and Brenner& Condiff equations greatly underpredict the measured data for the EG-TNT nanofluids.The high-shear viscosity of EG-TNT nanofluids is much higher than that of the EG-TiO2nanofluids containing spherical nanoparticles,indicating a strong particle shape effect on the shear viscosity of nanofluids.Although the shear-thinning behaviour of the nanofluids could be partially attributed to the structuring of interacting rod-like particles,the large deviation between the measured high-shear viscosity and the predicted ones by the Brenner&Condiff equation cannot fully be interpreted.In the following,an attempt is made to explain the experimental observations from the viewpoint of aggregation of nanaoparticles,which have been shown to play a key role in thermal behaviour of nanofluids in recent studies(Wang et al. 2003;Xuan et al.2003;Nan et al.1997;Prasher et al. 2006a,b;Keblinski et al.2005).Such an approach is also supported by the SEM and dynamic lightTable1Empirical constants for Eq.(1)a Maximum discrepancies;b Minimum discrepancies Concentration(wt%)A B C MaxD a(%)MinD b(%)0.0-3.21140.86973-154.570.62-1.440.5-3.42790.94425-148.490.93-0.471.0-2.94780.81435-159.14 1.11-0.692.0-2.2930.65293-174.57 1.64-0.694.0-2.63750.7574-165.820.99-0.948.0-2.73140.93156-145.010.88-1.57scattering analyses,which,as mentioned before, show clear evidence of particle aggregation.According to the modified Krieger-Dougherty equation(Goodwin and Hughes2000;Wang et al. 2003;Xuan et al.2003;Nan et al.1997),the relative viscosity of nanofluids,g r,is given as:g r¼1Àu a=u mðÞÀ½g u mð4Þwhere u m is the maximum concentration at which the flow can occur and u a is the effective volume fraction of aggregates given by u a¼u=u ma with u ma the maximum packing fraction of aggregates.As aggre-gates do not have constant packing throughout the structure,the packing density is assumed to change with radial position according to the power law with a constant index(D).As a result,u a is given as u a¼uða a=aÞ3ÀD;with a a and a,the effective radii of aggregates and primary nanoparticles,respectively. The term D is also referred as the fractal index meaning the extent of changes in the packing fraction from the centre to the edge of the aggregates.Typical values of D are given in normal textbook as D= 1.8–2.5for diffusion limited aggregation(DLA)and D=2.0–2.2for reaction limited aggregation(RLA); see for example Goodwin and Hughes(2000).For nanofluids containing spherical nanoparticles,the value of D has been shown experimentally and numerically to be between1.6and1.8(Wang et al. 2003,Xuan et al.2003)and between1.8and2.3, respectively(Waite et al.2001).A typical value of 1.8is suggested for nanofluids made of spherical nanoparticles(Prasher et al.2006a,b).However,little research has been found on the fractal index for nanofluids containing rod-like nanoparticles.The colloid science literature suggests a fractal index of 1.5–2.45for colloidal suspensions depending on the type of aggregation,chemistry environment,particle size and shape and shearflow conditions(Haas et al. 1993;Mohraz et al.2004;Hobbie and Fry2006; Micali et al.2006;Lin et al.2007).In a recent study, Mohraz et al.(2004)investigated the effect of monomer geometry on the fractal structure of colloi-dal rod aggregates.They found that the fractal index is a non-linear function of the monomer aspect ratio with the D increasing from*1.80to*2.3when the aspect ratio of the rod-like nanoparticles increases from1.0to30.6.Based on the above,a value of D=2.1is taken for nanofluids used in this work (Mohraz et al.2004,Lin et al.2007).Although the fractal model may appear to simplify the complexity of microstructures in aggregating systems containing rod-like particles,excellent agreement between the model prediction and experimental measurements exists when a a/a=9.46;see Fig.4.Here the aggregates are assumed to formflow units of an ellipsoidal shape with an effective aspect ratio of r a¼L a=b a;where L a and b a are the effective length and diameter,respectively.In the calculation,a typical value of u m of0.3is taken(Barnes et al.1989),and the intrinsic viscosity[g]is calculated by Eq.(3).It is to be noted that the aggregate size thatfits well to the rheological data(Fig.4)is consistent with the particle size analyses using both the SEM and the Malvern Nanosizer.A comparison between the EG-TNT data (a a/a=9.46,D=2.1,u m=0.30)and the EG-TiO2 data(a a/a=3.34,D=1.8,u m=0.605)(Chen et al. 2007b)in Fig.4suggests that the larger aggregate size in TNT nanofluids be an important factor responsible for the stronger shear-thinning behaviour and higher shear viscosity of TNT nanofluids.An inspection of Eq.(4)indicates that the effec-tive volume fraction u a u a¼u a a=aðÞ3ÀDis much higher than the actual volume fraction(u).This leads to the experimentally observed high-shear viscosity even for very dilute nanofluids,according to the classification discussed before.As a consequence,the demarcations defining the dilute and semi-concen-trated dispersions should be changed by using the effective volume fraction.The model discussed above can also provide a macroscopic explanation for the temperature indepen-dence of the high-shear viscosity.From Eq.(4),one can see that the relative high-shear viscosity depends on three parameters,the maximum volume fraction, u m,the effective volume fraction,u a and the intrinsic viscosity,[g].For a given nanofluid at a temperature not far from the ambient temperature,the three parameters are independent of temperature and hence the little temperature dependence of the relative shear viscosity.Microscopically,as explained before,the temperature-independent behaviour is due to negligi-ble Brownian diffusion compared with convection in high-shearflows.To further illustrate if the proposed aggregation mechanism is adequate,it is used to predict the effective thermal conductivity of the nanofluids by using the following conventional Hamilton–Crosser model(H–C model)(Hamilton and Crosser1962):k=k0¼k pþðnÀ1Þk0ÀðnÀ1Þuðk0Àk pÞk pþðnÀ1Þk0þuðk0Àk pÞð5Þwhere k and k0are,respectively,the thermal conductivities of nanofluids and base liquid,n is the shape factor given by n=3/w with w the surface area based sphericity.For TNT used in this work,the sphericity w is estimated as0.6(Hamilton and Crosser1962).For suspensions of aggregates,the above equation takes the following form:k=k0¼k aþðnÀ1Þk0ÀðnÀ1Þu aðk0Àk aÞa0a0að6Þwhere k a is the thermal conductivity of aggregates.To calculate k a,Eq.(6)is combined with the following Nan’s model(Nan et al.2003)for randomly dispersed nanotube-based composites:k a=k0¼3þu in½2b xð1ÀL xÞþb zð1ÀL zÞ3Àu in½2b x L xþb z L zð7Þwhere/in is the solid volume fraction of aggregates, b x¼ðk xÀk0Þ=½k mþL xðk tÀk mÞ and b z¼ðk zÀk0Þ=½k mþL zðk tÀk mÞ with k x,k m and k t being the thermal conductivities of nanotubes along trans-verse and longitudinal directions and isotropic thermal conductivity of the nanotube,respectively. In this work,k x,k m and k t are taken the same value as k p for afirst order of approximation due to lack of experimental data,and L x and L z are geometrical factors dependent on the nanotube aspect ratio given by L x¼0:5r2=ðr2À1ÞÀ0:5r coshÀ1r=ðr2À1Þ3=2 and L z¼1À2L x:Figure5shows the experimental results together with predictions by the original H–C model(Eq.5) and revised H–C model(Eq.6).Here the experiment data were obtained using a KD2thermal property meter(Labcell,UK)(Murshed et al.2005;Chen et al. 2008).One can see that the measured thermal conductivity is much higher than the prediction by the conventional H–C model(Eq.5),whereas the modified H–C model taking into account the effect of aggregation(Eq.6)agrees very well with the exper-imental data.The above results suggest that nanoparticle aggregates play a key role in the enhancement of thermal conductivity of nanofluids. The results also suggest that one could use the rheology data,which contain information of particle structuring in suspensions,for the effective thermal conductivity prediction.In summary,we have shown that EG-TNT nano-fluids are non-Newtonian exhibiting shear-thinning behaviour over20–60°C and a particle mass concen-tration range of0–8%,in contrast to the Newtonian behaviour for EG-TiO2nanofluids containing spher-ical particles.The non-Newtonian shear-thinning behaviour becomes stronger at higher temperatures or higher concentrations.For a given particle concen-tration,there exists a certain shear rate(e.g.*10s-1 for4wt%)below which the viscosity increases with increasing temperature,whereas the reverse occurs above such a shear rate.The normalised high-shearviscosity with respect to the base liquid viscosity, however,is found to be independent of temperature. These observations have not been reported in the literature for nanofluids.Further analyses suggest that the temperature effects are due to the shear-depen-dence of the relative contributions to the viscosity of the Brownian diffusion and convection.The analyses also suggest that a combination of particle aggregation and particle shape effects is the mechanism for the observed high-shear rheological behaviour,which is supported not only by the particle size measurements but also by the thermal conductivity measurements and analyses using a combination of the H–C and Nan’s models.The results of this work also indicate that one could use the information of aggregation from the rheological experiments for predicting the effec-tive thermal conductivity of nanofluids. Acknowledgement The work was partially supported by UK EPSRC under grants EP/E00041X/1and EP/F015380/1.ReferencesAbdulagatov MI,Azizov ND(2006)Experimental study of the effect of temperature,pressure and concentration on the viscosity of aqueous NaBr solutions.J Solut Chem 35(5):705–738.doi:10.1007/s10953-006-9020-6Barnes HA,Hutton JF,Walters K(1989)An introduction to rheology.Elsevier Science B.V.,NetherlandsBavykin DV,Parmon VN,Lapkin AA,Walsh FC(2004)The effect of hydrothermal conditions on the mesoporous struc-ture of TiO2nanotubes.J Mater Chem14(22):3370–3377 Bird RB,Steward WE and Lightfoot EN(2002)Transport Phenomena,2nd edn.Wiley,New YorkBrenner H,Condiff DW(1974)Transport mechanics in sys-tems of orientable particles,Part IV.Convective Transprort.J Colloid Inter Sci47(1):199–264Chen HS,Ding YL,He YR,Tan CQ(2007a)Rheological behaviour of ethylene glycol based titania nanofluids.Chem Phys Lett444(4–6):333–337Chen HS,Ding YL,Tan CQ(2007b)Rheological behaviour of nanofluids.New J Phys9(367):1–25Chen HS,Yang W,He YR,Ding YL,Zhang LL,Tan CQ,Lapkin AA,Bavykin DV(2008)Heat transfer andflow behaviour of aqueous suspensions of titanate nanotubes under the laminar flow conditions.Powder Technol183:63–72Choi SUS(1995)Enhancing thermal conductivity offluids with nanoparticles In:Siginer DA,Wang HP(eds) Developments applications of non-newtonianflows,FED-vol231/MD-vol66.ASME,New York,pp99–105 Chow TS(1993)Viscosities of concentrated dispersions.Phys Rev E48:1977–1983Das SK,Putra N,Roetzel W(2003)Pool boiling characteristics of nano-fluids.Int J Heat Mass Transfer46:851–862Ding YL,Alias H,Wen DS,Williams RA(2006)Heat transfer of aqueous suspensions of carbon nanotubes(CNT nanofluids).Int J Heat Mass Transf49(1–2):240–250 Doi M,Edwards SF(1978a)Dynamics of rod-like macro-molecules in concentrated solution,Part1.J Colloid Sci 74:560–570Doi M,Edwards SF(1978b)Dynamics of rod-like macro-molecules in concentrated solution,Part2.J Colloid Sci 74:918–932Eastman JA,Choi SUS,Li S,Yu W,Thompson LJ(2001) Anomalously increased effective thermal conductivities of ethylene glycol-based nanofluids containing copper nanoparticles.Appl Phys Lett78:718–720Egres RG,Wagner NJ(2005)The rheology and microstructure of acicular precipitated calcium carbonate colloidal sus-pensions through the shear thickening transition.J Rheol 49:719–746Einstein A(1906)Eine neue Bestimmung der Molekul-dimension(a new determination of the molecular dimensions).Annal der Phys19(2):289–306Einstein A(1911)Berichtigung zu meiner Arbeit:Eine neue Bestimmung der Molekul-dimension(correction of my work:a new determination of the molecular dimensions).Ann der Phys34(3):591–592Goodwin JW,Hughes RW(2000)Rheology for chemists—an introduction.The Royal Society of Chemistry,UK Haas W,Zrinyi M,Kilian HG,Heise B(1993)Structural analysis of anisometric colloidal iron(III)-hydroxide par-ticles and particle-aggregates incorporated in poly(vinyl-acetate)networks.Colloid Polym Sci271:1024–1034 Hamilton RL,Crosser OK(1962)Thermal Conductivity of heterogeneous two-component systems.I&EC Fundam 125(3):187–191He YR,Jin Y,Chen HS,Ding YL,Cang DQ(2007)Heat transfer andflow behaviour of aqueous suspensions of TiO2nanoparticles(nanofluids)flowing upward through a vertical pipe.Int J Heat Mass Transf50(11–12):2272–2281Hobbie EK,Fry DJ(2006)Nonequilibrium phase diagram of sticky nanotube suspensions.Phys Rev Lett97:036101 Keblinski P,Eastman JA and Cahill DG(2005)Nanofluids for thermal transport,Mater Today,June Issue,36–44 Krishnamurthy S,Lhattacharya P,Phelan PE,Prasher RS (2006)Enhanced mass transport of in nanofluids.Nano Lett6(3):419–423Kwak K,Kim C(2005)Viscosity and thermal conductivity of copper oxide nanofluid dispersed in ethylene glycol.Korea-Aust Rheol J17(2):35–40Larson RG(1999)The structure and rheology of complex fluids.Oxford University Press,New YorkLarson RG(2005)The rheology of dilute solutions offlexible polymers:progress and problems.J Rheol49:1–70Lee D,Kim J,Kim B(2006)A new parameter to control heat transport in nanofluids:surface charge state of the particle in suspension.J Phys Chem B110:4323–4328Lin JM,Lin TL,Jeng U,Zhong Y,Yeh C,Chen T(2007) Fractal aggregates of fractal aggregates of the Pt nano-particles synthesized by the polyol process and poly(N-vinyl-2-pyrrolidone)reduction.J Appl Crystallogr40: s540–s543Micali N,Villari V,Castriciano MA,Romeo A,Scolaro LM (2006)From fractal to nanorod porphyrin J-aggregates.Concentration-induced tuning of the aggregate size.J Phys Chem B110:8289–8295Mohraz A,Moler DB,Ziff RM,Solomon MJ(2004)Effect of monomer geometry on the fractal structure of colloidal rod aggregates.Phys Rev Lett92:155503Murshed SMS,Leong KC,Yang C(2005)Enhanced thermal conductivity of TiO2-water based nanofluids.Int J Therm Sci44:367–373Nan CW,Birringer R,Clarke DR,Gleiter H(1997)Effective thermal conductivity of particulate composites with inter-facial thermal resistance.J Appl Phys81(10):6692–6699 Nan CW,Shi Z,Lin Y(2003)A simple model for thermal conductivity of carbon nanotube-based composites.Chem Phys Lett375(5–6):666–669Pak BC,Cho YI(1998)Hydrodynamic and heat transfer study of dispersedfluids with submicron metallic oxide parti-cles.Exp Heat Transf11:150–170Prasher R,Song D,Wang J(2006a)Measurements of nanofluid viscosity and its implications for thermal applications.Appl Phys Lett89:133108-1-3Prasher R,Phelan PE,Bhattacharya P(2006b)Effect of aggregation kinetics on thermal conductivity of nanoscale colloidal solutions(nanofluids).Nano Lett6(7):1529–1534Russel WB,Saville DA,Scholwater WR(1991)Colloidal dispersions.Cambridge University Press,Cambridge Waite TD,Cleaver JK,Beattie JK(2001)Aggregation kinetics and fractal structure of gamma-alumina assemblages.J Colloid Interface Sci241:333–339Wang XQ,Mujumdar AS(2007)Heat transfer characteristics of nanofluids:a review.Int J Therm Sci46:1–19Wang XW,Xu XF,Choi SUS(1999)Thermal conductivity of nanoparticle–fluid mixture.J Thermphys Heat Transf 13:474Wang BX,Zhou LP,Peng XF(2003)A fractal model for predicting the effective thermal conductivity of liquid with suspension of nanoparticles.Int J Heat Mass Transf 46:2665–2672Wasan DT,Nikolov AD(2003)Spreading of nanofluids on solids.Nature423(6936):156–159Wen DS,Ding YL(2005)Formulation of nanofluids for natural convective heat transfer applications.Int J Heat Fluid Flow26:855–864Xuan YM,Li Q,Hu J,WF(2003)Aggregation structure and thermal conductivity of nanofluids.AIChE J49(4):1038–1043Yang Y,Zhong ZG,Grulke EA,Anderson WB,Wu G(2005) Heat transfer properties of nanoparticle-in-fluid dispersion (nanofluids)in laminarflow.Int J Heat Mass Transf 48:1107–1116。

International Journal of Computer Vision 65(1/2), 43–72, 2005c 2005Springer Science +Business Media, Inc. Manufactured in The Netherlands.DOI:10.1007/s11263-005-3848-xA Comparison of Affine Region DetectorsK.MIKOLAJCZYKUniversity of Oxford,OX13PJ,Oxford,United Kingdomkm@T.TUYTELAARSUniversity of Leuven,Kasteelpark Arenberg10,3001Leuven,Belgiumtuytelaa@esat.kuleuven.beC.SCHMIDINRIA,GRAVIR-CNRS,655,av.de l’Europe,38330,Montbonnot,Franceschmid@inrialpes.frA.ZISSERMANUniversity of Oxford,OX13PJ,Oxford,United Kingdomaz@J.MATASCzech Technical University,Karlovo Namesti13,12135,Prague,Czech Republicmatas@cmp.felk.cvut.czF.SCHAFFALITZKY AND T.KADIRUniversity of Oxford,OX13PJ,Oxford,United Kingdomfsm@tk@L.V AN GOOLUniversity of Leuven,Kasteelpark Arenberg10,3001Leuven,Belgiumvangool@esat.kuleuven.beReceived August20,2004;Revised May3,2005;Accepted May11,2005First online version published in January,2006Abstract.The paper gives a snapshot of the state of the art in affine covariant region detectors,and compares their performance on a set of test images under varying imaging conditions.Six types of detectors are included: detectors based on affine normalization around Harris(Mikolajczyk and Schmid,2002;Schaffalitzky and Zisserman,2002)and Hessian points(Mikolajczyk and Schmid,2002),a detector of‘maximally stable extremal regions’,proposed by Matas et al.(2002);an edge-based region detector(Tuytelaars and Van Gool,1999) and a detector based on intensity extrema(Tuytelaars and Van Gool,2000),and a detector of‘salient regions’,44Mikolajczyk et al.proposed by Kadir,Zisserman and Brady(2004).The performance is measured against changes in viewpoint,scale, illumination,defocus and image compression.The objective of this paper is also to establish a reference test set of images and performance software,so that future detectors can be evaluated in the same framework.Keywords:affine region detectors,invariant image description,local features,performance evaluation1.IntroductionDetecting regions covariant with a class of transforma-tions has now reached some maturity in the computer vision literature.These regions have been used in quite varied applications including:wide baseline matching for stereo pairs(Baumberg,2000;Matas et al.,2002;Pritchett and Zisserman,1998;Tuytelaars and Van Gool,2000),reconstructing cameras for sets of disparate views(Schaffalitzky and Zisserman, 2002),image retrieval from large databases(Schmid and Mohr,1997;Tuytelaars and Van Gool,1999), model based recognition(Ferrari et al.,2004;Lowe, 1999;Obdrˇz´a lek and Matas,2002;Rothganger et al., 2003),object retrieval in video(Sivic and Zisserman, 2003;Sivic et al.,2004),visual data mining(Sivic and Zisserman,2004),texture recognition(Lazebnik et al.,2003a,b),shot location(Schaffalitzky and Zisserman,2003),robot localization(Se et al.,2002) and servoing(Tuytelaars et al.,1999),building panoramas(Brown and Lowe,2003),symmetry detection(Turina et al.,2001),and object categoriza-tion(Csurka et al.,2004;Dorko and Schmid,2003; Fergus et al.,2003;Opelt et al.,2004).The requirement for these regions is that they should correspond to the same pre-image for dif-ferent viewpoints,i.e.,their shape is notfixed but automatically adapts,based on the underlying image intensities,so that they are the projection of the same 3D surface patch.In particular,consider images from two viewpoints and the geometric transformation between the images induced by the viewpoint change. Regions detected after the viewpoint change should be the same,modulo noise,as the transformed versions of the regions detected in the original image–image transformation and region detection commute.As such,even though they have often been called invariant regions in the literature(e.g.,Dorko and Schmid,2003;Lazebnik et al.,2003a;Sivic and Zisserman,2004;Tuytelaars and Van Gool,1999),in principle they should be termed covariant regions since they change covariantly with the transformation. The confusion probably arises from the fact that,even though the regions themselves are covariant,the nor-malized image pattern they cover and the feature de-scriptors derived from them are typically invariant. Note,our use of the term‘region’simply refers to a set of pixels,i.e.any subset of the image.This differs from classical segmentation since the region bound-aries do not have to correspond to changes in image appearance such as colour or texture.All the detectors presented here produce simply connected regions,but in general this need not be the case.For viewpoint changes,the transformation of most interest is an affinity.This is illustrated in Fig.1. Clearly,a region withfixed shape(a circular exam-ple is shown in Fig.1(a)and(b))cannot cope with the geometric deformations caused by the change in view-point.We can observe that the circle does not cover the same image content,i.e.,the same physical surface. Instead,the shape of the region has to be adaptive,or covariant with respect to affinities(Fig.1(c)–close-ups shown in Fig.1(d)–(f)).Indeed,an affinity is suffi-cient to locally model image distortions arising from viewpoint changes,provided that(1)the scene sur-face can be locally approximated by a plane or in case of a rotating camera,and(2)perspective effects are ignored,which are typically small on a local scale any-way.Aside from the geometric deformations,also pho-tometric deformations need to be taken into account. These can be modeled by a linear transformation of the intensities.To further illustrate these issues,and how affine covariant regions can be exploited to cope with the geometric and photometric deformation between wide baseline images,consider the example shown in Fig.2.Unlike the example of Fig.1(where a circular region was chosen for one viewpoint)the elliptical image regions here are detected independently in each viewpoint.As is evident,the pre-images of these affineA Comparison of Affine Region Detectors45Figure 1.Class of transformations needed to cope with viewpoint changes.(a)First viewpoint;(b,c)second viewpoint.Fixed size circular patches (a,b)clearly do not suffice to deal with general viewpoint changes.What is needed is an anisotropic rescaling,i.e.,an affinity (c).Bottom row shows close-up of the images of the toprow.Figure 2.Affine covariant regions offer a solution to viewpoint and illumination changes.First row:one viewpoint;second row:other viewpoint.(a)Original images,(b)detected affine covariant regions,(c)close-up of the detected regions.(d)Geometric normalization to circles.The regions are the same up to rotation.(e)Photometric and geometric normalization.The slight residual difference in rotation is due to an estimation error.covariant regions correspond to the same surface region.Given such an affine covariant region,it is then possible to normalize against the geometric and photometric deformations (shown in Fig.2(d),(e))and to obtain a viewpoint and illumination invariant description of the intensity pattern within the region.In a typical matching application,the regions are used as follows.First,a set of covariant regions is46Mikolajczyk et al.detected in an image.Often a large number,perhaps hundreds or thousands,of possibly overlapping regions are obtained.A vector descriptor is then asso-ciated with each region,computed from the intensity pattern within the region.This descriptor is chosen to be invariant to viewpoint changes and,to some extent, illumination changes,and to discriminate between the regions.Correspondences may then be established with another image of the same scene,byfirst detect-ing and representing regions(independently)in the new image;and then matching the regions based on their descriptors.By design the regions commute with viewpoint change,so by design,corresponding regions in the two images will have similar(ideally identical) vector descriptors.The benefits are that correspon-dences can then be easily established and,since there are multiple regions,the method is robust to partial occlusions.This paper gives a snapshot of the state of the art in affine covariant region detection.We will describe and compare six methods of detecting these regions on images.These detectors have been designed and implemented by a number of researchers and the comparison is carried out using binaries supplied by the authors.The detectors are:(i)the ‘Harris-Affine’detector(Mikolajczyk and Schmid, 2002,2004;Schaffalitzky and Zisserman,2002); (ii)the‘Hessian-Affine’detector(Mikolajczyk and Schmid,2002,2004);(iii)the‘maximally stable extremal region’detector(or MSER,for short)(Matas et al.,2002,2004);(iv)an edge-based region detec-tor(Tuytelaars and Van Gool,1999,2004)(referred to as EBR);(v)an intensity extrema-based region detector(Tuytelaars and Van Gool,2000,2004) (referred to as IBR);and(vi)an entropy-based region detector(Kadir et al.,2004)(referred to as salient regions).To limit the scope of the paper we have not included methods for detecting regions which are covariant only to similarity transformations(i.e.,in particular scale), such as(Lowe,1999,2004;Mikolajczyk and Schmid, 2001;Mikolajczyk et al.,2003),or other methods of computing affine invariant descriptors,such as image lines connecting interest points(Matas et al.,2000; Tell and Carlson,2000,2002),or invariant vertical line segments(Goedeme et al.,2004).Also the detectors proposed by Lindeberg and G˚a rding(1997)and Baum-berg(2000)have not been included,as they come very close to the Harris-Affine and Hessian-Affine detectors.The six detectors are described in Section2.They are compared on the data set shown in Fig.9.This data set includes structured and textured scenes as well as different types of transformations:viewpoint changes,scale changes,illumination changes,blur and JPEG compression.It is described in more detail in Section3.Two types of comparisons are carried out. First,in Section10,the repeatability of the detector is measured:how well does the detector determine cor-responding scene regions?This is measured by com-paring the overlap between the ground truth and de-tected regions,in a manner similar to the evaluation test used in Mikolajczyk and Schmid(2002),but with special attention paid to the effect of the different scales (region sizes)of the various detectors’output.Here, we also measure the accuracy of the regions’shape, scale and localization.Second,the distinctiveness of the detected regions is assessed:how distinguishable are the regions detected?Following(Mikolajczyk and Schmid,2003,2005),we use the SIFT descriptor de-veloped by Lowe(1999),which is an128-dimensional vector,to describe the intensity pattern within the im-age regions.This descriptor has been demonstrated to be superior to others used in literature on a number of measures(Mikolajczyk and Schmid,2003).Our intention is that the images and tests de-scribed here will be a benchmark against which fu-ture affine covariant region detectors can be assessed. The images,Matlab code to carry out the performance tests,and binaries of the detectors are available from /∼vgg/research/affine. 2.Affine Covariant DetectorsIn this section we give a brief description of the six re-gion detectors used in the comparison.Section2.1de-scribes the related methods Harris-Affine and Hessian-Affine.Sections2.2and2.3describe methods for detecting edge-based regions and intensity extrema-based regions.Finally,Sections2.4and2.5describe MSER and salient regions.For the purpose of the comparisons the output re-gion of all detector types are represented by a common shape,which is an ellipse.Figures3and4show the el-lipses for all detectors on one pair of images.In order not to overload the images,only some of the corre-sponding regions that were actually detected in both images have been shown.This selection is obtained by increasing the threshold.A Comparison of Affine Region Detectors47Figure3.Regions generated by different detectors on corresponding sub-parts of thefirst and third graffiti images of Fig.9(a).The ellipses show the original detection size.In fact,for most of the detectors the output shape is an ellipse.However,for two of the de-tectors(edge-based regions and MSER)it is not, and information is lost by this representation,as ellipses can only be matched up to a rotational degree of freedom.Examples of the original re-gions detected by these two methods are given in Fig.5.These are parallelogram-shaped regions for the edge-based region detector,and arbitrarily shaped regions for the MSER detector.In the following the representing ellipse is chosen to have the same first and second moments as the originally detected region,which is an affine covariant construction method.48Mikolajczyk etal.Figure 4.Regions generated by different detectors continued.2.1.Detectors Based on Affine Normalization—Harris-Affine &Hessian-AffineWe describe here two related methods which detect interest points in scale-space,and then determine an elliptical region for each point.Interest points are either detected with the Harris detector or with a detector based on the Hessian matrix.In both casesscale-selection is based on the Laplacian,and the shape of the elliptical region is determined with the second moment matrix of the intensity gradient (Baumberg,2000;Lindeberg and G˚a rding,1997).The second moment matrix,also called the auto-correlation matrix,is often used for feature detection or for describing local image structures.Here it is used both in the Harris detector and the ellipticalA Comparison of Affine Region Detectors49Figure5.Originally detected region shapes for the regions shown in Figs.3(c)and4(b). shape estimation.This matrix describes the gradientdistribution in a local neighbourhood of a point:M=µ(x,σI,σD)= µ11µ12µ21µ22=σ2D g(σI)∗I2x(x,σD)I x I y(x,σD)I x I y(x,σD)I2y(x,σD)(1)The local image derivatives are computed with Gaussian kernels of scaleσD(differentiation scale). The derivatives are then averaged in the neighbourhood of the point by smoothing with a Gaussian window of scaleσI(integration scale).The eigenvalues of this matrix represent two principal signal changes in a neighbourhood of the point.This property enables the extraction of points,for which both curvatures are significant,that is the signal change is significant in orthogonal directions.Such points are stable in arbitrary lighting conditions and are representative of an image.One of the most reliable interest point detectors,the Harris detector(Harris and Stephens, 1988),is based on this principle.A similar idea is explored in the detector based on the Hessian matrix:H=H(x,σD)=h11h12h21h22=I xx(x,σD)I xy(x,σD)I xy(x,σD)I yy(x,σD)(2)The second derivatives,which are used in this matrix give strong responses on blobs and ridges.The regions are similar to those detected by a Laplacian operator (trace)(Lindeberg,1998;Lowe,1999)but a function based on the determinant of the Hessian matrix penal-izes very long structures for which the second deriva-tive in one particular orientation is very small.A local maximum of the determinant indicates the presence of a blob structure.To deal with scale changes a scale selection method(Lindeberg,1998)is applied.The idea is to select the characteristic scale of a local structure, for which a given function attains an extremum over scales(see Fig.6).The selected scale is characteristic in the quantitative sense,since it measures the scale50Mikolajczyk etal.Figure 6.Example of characteristic scales.Top row shows images taken with different zoom.Bottom row shows the responses of the Laplacian over scales.The characteristic scales are 10.1and 3.9for the left and right image,respectively.The ratio of scales corresponds to the scale factor (2.5)between the two images.The radius of displayed regions in the top row is equal to 3times the selected scales.at which there is maximum similarity between the feature detection operator and the local image struc-tures.The size of the region is therefore selected in-dependently of image resolution for each point.The Laplacian operator is used for scale selection in both detectors since it gave the best results in the ex-perimental comparison in Mikolajczyk and Schmid (2001).Given the set of initial points extracted at their char-acteristic scales we can apply the iterative estimation of elliptical affine region (Lindeberg and G˚a rding,1997).The eigenvalues of the second moment matrix are used to measure the affine shape of the point neighbourhood.To determine the affine shape,we find the transforma-tion that projects the affine pattern to the one with equal eigenvalues.This transformation is given by the square root of the second moment matrix M 1/2.If the neigh-bourhood of points x R and x L are normalized by trans-formations x R =M 1/2R x R and x L =M 1/2L x L ,respec-tively,the normalized regions are related by a simple rotation x L =R x R (Baumberg,2000;Lindeberg andG˚a rding,1997).The matrices M Land M R computed in the normalized frames are equal to a rotation matrix (see Fig.7).Note that rotation preserves the eigen-value ratio for an image patch,therefore,the affine deformation can be determined up to a rotation fac-tor.This factor can be recovered by other methods,for example normalization based on the dominant gradi-ent orientation (Lowe,1999;Mikolajczyk and Schmid,2002).The estimation of affine shape can be applied to any initial point given that the determinant of the second moment matrix is larger than zero and the signal to noise ratio is insignificant for this point.We can there-fore use this technique to estimate the shape of initial regions provided by the Harris and Hessian based de-tector.The outline of the iterative region estimation:1.Detect initial region with Harris or Hessian detector and select the scale.2.Estimate the shape with the second moment matrix3.Normalize the affine region to the circular one4.Go to step 2if the eigenvalues of the second moment matrix for new point are not equal.Examples of Harris-Affine and Hessian-Affine re-gions are displayed on Fig.3(a)and (b).2.2.An Edge-Based Region DetectorWe describe here a method to detect affine covariant regions in an image by exploiting the edges present in the image.The rationale behind this is that edges are typically rather stable features,that can be detected over a range of viewpoints,scales and/or illumination changes.Moreover,by exploiting the edge geometry,the dimensionality of the problem can be significantly reduced.Indeed,as will be shown next,the 6D search problem over all possible affinities (or 4D,once theA Comparison of Affine Region Detectors51Figure 7.Diagram illustrating the affine normalization using the second moment matrices.Image coordinates are transformed with matrices M −1/2L and M −1/2R .center point is fixed)can further be reduced to a one-dimensional problem by exploiting the nearby edges geometry.In practice,we start from a Harris corner point p (Harris and Stephens,1988)and a nearby edge,extracted with the Canny edge detector (Canny,1986).To increase the robustness to scale changes,these basic features are extracted at multiple scales.Two points p 1and p 2move away from the corner in both directions along the edge,as shown in Fig.8(a).Their relative speed is coupled through the equality of relative affine invariant parameters l 1and l 2:l i =absp i (1)(s i )p −p i (s i ) ds i (3)with s i an arbitrary curve parameter (in both direc-tions),p i (1)(s i )the first derivative of p i (s i )with respectto s i ,abs()the absolute value and |...|the determi-nant.This condition prescribes that the areas between the joint p ,p 1 and the edge and between the joint p ,p 2 and the edge remain identical.This is an affine invariant criterion indeed.From now on,we simply use l when referring to l 1=l 2.For each value l ,the two points p 1(1)and p 2(1)together with the corner p define a parallelogram (l ):the parallelogram spanned by the vectors p 1(l )−p and p 2(l )−p .This yields a one dimensional family of parallelogram-shaped regions as a function of l .From this 1D family we select one (or a few)parallelogram for which the following photometric quantities of the texture go through an extremum.Inv 1=abs|p 1−p g p 2−p g ||p −p 1p −p 2| M 100 M 200M 000−(M 100)2Inv 2=abs|p −p g q −p g |1p −p 2 M 100 M 200M 000−(M 100)2Figure 8.Construction methods for EBR and IBR.(a)The edge-based region detector starts from a corner point p and exploits nearby edgeinformation;(b)The intensity extrema-based region detector starts from an intensity extremum and studies the intensity pattern along rays emanating from this point.52Mikolajczyk et al.with M n pq =I n (x ,y )x p y q dxdy(4)p g = M 110M 100,M 101M 100with M n pq the n th order ,(p +q )th degree momentcomputed over the region (l ),p g the center of gravity of the region,weighted with intensity I (x ,y ),and q the corner of the parallelogram opposite to the corner point p (see Fig.8(a)).The second factor in these formula has been added to ensure invariance under an intensity offset.In the case of straight edges,the method described above cannot be applied,since l =0along the entire edge.Since intersections of two straight edges occur quite often,we cannot simply neglect this case.To cir-cumvent this problem,the two photometric quantities given in Eq.(4)are combined and locations where both functions reach a minimum value are taken to fix the parameters s 1and s 2along the straight edges.Moreover,instead of relying on the correct detection of the Harris corner point,we can simply use the straight lines intersection point instead.A more detailed expla-nation of this method can be found in Tuytelaars and Van Gool (1999,2004).Examples of detected regions are displayed in Fig.5(b).For easy comparison in the context of this paper,the parallelograms representing the invariant regions are replaced by the enclosed ellipses,as shown in Fig.4(b).However,in this way the orientation-information is lost,so it should be avoided in a practical application,as discussed in the beginning of Section 2.2.3.Intensity Extrema-Based Region DetectorHere we describe a method to detect affine covariant regions that starts from intensity extrema (detected at multiple scales),and explores the image around them in a radial way,delineating regions of arbitrary shape,which are then replaced by ellipses.More precisely,given a local extremum in inten-sity,the intensity function along rays emanating from the extremum is studied,as shown in Fig.8(b).The following function is evaluated along each ray:f I (t )=abs (I (t )−I 0)max t 0abs (I (t )−I 0)dt t,d with t an arbitrary parameter along the ray,I (t )the intensity at position t ,I 0the intensity value at the ex-tremum and d a small number which has been added to prevent a division by zero.The point for which this function reaches an extremum is invariant under affine geometric and linear photometric transforma-tions (given the ray).Typically,a maximum is reached at positions where the intensity suddenly increases or decreases.The function f I (t )is in itself already in-variant.Nevertheless,we select the points where this function reaches an extremum to make a robust selec-tion.Next,all points corresponding to maxima of f I (t )along rays originating from the same local extremum are linked to enclose an affine covariant region (see Fig.8(b)).This often irregularly-shaped region is re-placed by an ellipse having the same shape moments up to the second order.This ellipse-fitting is again an affine covariant construction.Examples of detected re-gions are displayed in Fig.4(a).More details about this method can be found in Tuytelaars and Van Gool (2000,2004).2.4.Maximally Stable Extremal Region DetectorA Maximally Stable Extremal Region (MSER)is aconnected component of an appropriately thresholded image.The word ‘extremal’refers to the property that all pixels inside the MSER have either higher (bright extremal regions)or lower (dark extremal regions)intensity than all the pixels on its outer boundary.The ‘maximally stable’in MSER describes the property optimized in the threshold selection process.The set of extremal regions E ,i.e.,the set of all connected components obtained by thresholding,has a number of desirable properties.Firstly,a mono-tonic change of image intensities leaves E unchanged,since it depends only on the ordering of pixel intensi-ties which is preserved under monotonic transforma-tion.This ensures that common photometric changes modelled locally as linear or affine leave E unaffected,even if the camera is non-linear (gamma-corrected).Secondly,continuous geometric transformations pre-serve topology–pixels from a single connected compo-nent are transformed to a single connected component.Thus after a geometric change locally approximated by an affine transform,homography or even continuous non-linear warping,a matching extremal region will be in the transformed set E .Finally,there are no more extremal regions than there are pixels in the image.SoA Comparison of Affine Region Detectors53a set of regions was defined that is preserved under a broad class of geometric and photometric changes and yet has the same cardinality as e.g.the set offixed-sized square windows commonly used in narrow-baseline matching.Implementation Details.The enumeration of the set of extremal regions E is very efficient,almost linear in the number of image pixels.The enumeration pro-ceeds as follows.First,pixels are sorted by intensity. After sorting,pixels are marked in the image(either in decreasing or increasing order)and the list of growing and merging connected components and their areas is maintained using the union-find algorithm(Sedgewick, 1988).During the enumeration process,the area of each connected component as a function of intensity is stored.Among the extremal regions,the‘maximally stable’ones are those corresponding to thresholds were the relative area change as a function of relative change of threshold is at a local minimum.In other words, the MSER are the parts of the image where local bi-narization is stable over a large range of thresholds. The definition of MSER stability based on relative area change is only affine invariant(both photomet-rically and geometrically).Consequently,the process of MSER detection is affine covariant.Detection of MSER is related to thresholding,since every extremal region is a connected component of a thresholded image.However,no global or‘optimal’threshold is sought,all thresholds are tested and the stability of the connected components evaluated.The output of the MSER detector is not a binarized image. For some parts of the image,multiple stable thresholds exist and a system of nested subsets is output in this case.Finally we remark that the different sets of extremal regions can be defined just by changing the ordering function.The MSER described in this section and used in the experiments should be more precisely called intensity induced MSERs.2.5.Salient Region DetectorThis detector is based on the pdf of intensity values computed over an elliptical region.Detection proceeds in two steps:first,at each pixel the entropy of the pdf is evaluated over the three parameter family of el-lipses centred on that pixel.The set of entropy extrema over scale and the corresponding ellipse parameters are recorded.These are candidate salient regions.Second,the candidate salient regions over the entire image are ranked using the magnitude of the derivative of the pdf with respect to scale.The top P ranked regions are retained.In more detail,the elliptical region E centred on a pixel x is parameterized by its scale s(which specifies the major axis),its orientationθ(of the major axis), and the ratio of major to minor axesλ.The pdf of intensities p(I)is computed over E.The entropy H is then given byH=−Ip(I)log p(I)The set of extrema over scale in H is computed for the parameters s,θ,λfor each pixel of the image.For each extrema the derivative of the pdf p(I;s,θ,λ)with s is computed asW=s22s−1I∂p(I;s,θ,λ)∂s,and the saliency Y of the elliptical region is com-puted as Y=HW.The regions are ranked by their saliency Y.Examples of detected regions are displayed in Fig.4(c).More details about this method can be found in Kadir et al.(2004).3.The Image Data SetFigure9shows examples from the image sets used to evaluate the detectors.Five different changes in imag-ing conditions are evaluated:viewpoint changes(a) &(b);scale changes(c)&(d);image blur(e)&(f); JPEG compression(g);and illumination(h).In the cases of viewpoint change,scale change and blur,the same change in imaging conditions is applied to two different scene types.This means that the effect of changing the image conditions can be separated from the effect of changing the scene type.One scene type contains homogeneous regions with distinctive edge boundaries(e.g.graffiti,buildings),and the other con-tains repeated textures of different forms.These will be referred to as structured versus textured scenes re-spectively.In the viewpoint change test the camera varies from a fronto-parallel view to one with significant fore-shortening at approximately60degrees to the camera. The scale change and blur sequences are acquired by varying the camera zoom and focus respectively.。

Geometric ModelingGeometric modeling is a branch of mathematics that deals with the representation of objects in space. It is a fundamental tool in computer graphics, computer-aided design (CAD), and other applications that require the creation of 3D models. Geometric modeling involves the use of mathematical equations and algorithms to create and manipulate objects in space. In this essay, we will explore the different aspects of geometric modeling, including its history, applications, and challenges.The history of geometric modeling can be traced back to the early 19th century when mathematicians began to study the properties of curves and surfaces. In the early 20th century, the development of calculus and differential geometry led to the creation of new methods for representing complex objects in space. The introduction of computers in the mid-20th century revolutionized the field of geometric modeling, making it possible to create and manipulate 3D models with greater precision and ease.Today, geometric modeling is used in a wide range of applications, including computer graphics, animation, video games, virtual reality, and CAD. In computer graphics and animation, geometric modeling is used to create realistic 3D models of objects, characters, and environments. In video games, geometric modeling is used to create the game world and characters. In virtual reality, geometric modeling is used to create immersive environments that simulate real-world experiences. In CAD, geometric modeling is used to create precise 3D models of mechanical parts and assemblies.One of the biggest challenges in geometric modeling is the representation of complex shapes and surfaces. Many real-world objects, such as cars, airplanes, and human bodies, have complex shapes that are difficult to represent using simple geometric primitives such as spheres, cylinders, and cones. To overcome this challenge, researchers have developed advanced techniques such as NURBS (non-uniform rational B-splines), which allow for the creation of complex curves and surfaces by combining simple geometric primitives.Another challenge in geometric modeling is the optimization of models for efficient rendering and simulation. As the complexity of models increases, so doesthe computational cost of rendering and simulating them. To address this challenge, researchers have developed techniques such as level-of-detail (LOD) modeling,which involves creating multiple versions of a model at different levels of detail to optimize rendering and simulation performance.In conclusion, geometric modeling is a fundamental tool in computer graphics, animation, video games, virtual reality, and CAD. Its history can be traced backto the early 19th century, and its development has been closely tied to the advancement of mathematics and computer technology. Despite its many applications and successes, geometric modeling still faces challenges in the representation of complex shapes and surfaces, as well as the optimization of models for efficient rendering and simulation. As technology continues to advance, it is likely that new techniques and approaches will emerge to overcome these challenges and pushthe field of geometric modeling forward.。

Book one Unit1Read the following passage carefully and choose the best word or phrase given below to fill in each blank.Change the form where necessary.A recent survey of women in 20 large and medium-sized cities across the country revealed that about half of the respondents were happy with their marriages and relationships, while nearly 30 percent said they were bored (1)and 3.4 percent they were in agony.近日，一项针对全国20个大中城市女性的调查显示，约半数受访者对自己的婚姻和感情生活感到满意，近30%的人对此感到厌倦，3.4%的人感到痛苦。

3 percent said they were worried about their relationships and 12 percent said the did not know how to describe their mixed feelings.The Huakun Women Survey Center, a(n) affiliate (2)of the All-China Women’s Federation, conduct ed (3)the survey of 2,000 women aged between 20 and 60 at the end of last year. Altogether 1,955 valid (4)questionnaires were collected. 3%的人说他们担心自己的恋情，12%的人说他们不知道如何描述自己复杂的感情。

3D shape measurement technique formultiple rapidly moving objectsYajun Wang,Song Zhang,∗and James H.OliverDepartment of Mechanical Engineering,Iowa State University,Ames,Iowa50011,USA∗song@Abstract:Recently proposed binary defocusing techniques have led toultrafast speed3D shape measurement,but they are generally limited tomeasurement of a single object at a time.Introducing additional gray codedpatterns for point-by-point phase unwrapping could permit simultaneousmultiple-object measurement.However,when the objects are movingrapidly,the displacement between thefirst captured pattern and the last canbe significant,and pose challenges related to the precisely designed graycodes.This paper presents a new phase unwrapping strategy that combinesthe conventional spatial phase unwrapping with the gray code to resolvemotion related phase unwrapping problems.A system with a speed of5,000Hz was developed to verify the performance of the proposed technique.©2011Optical Society of AmericaOCIS codes:(320.7100)Ultrafast measurements;(110.6880)Three-dimensional image acqui-sition;(120.5050)Phase measurement.References and links1.S.Gorthi and P.Rastogi,“Fringe projection techniques:whither we are?”ser.Eng.48,133–140(2010).2.S.Zhang,“Recent progresses on real-time3-D shape measurement using digital fringe projection techniques,”ser Eng.48(2),149–158(2010).3.S.Lei and S.Zhang,“Flexible3-D shape measurement using projector defocusing,”Opt.Lett.34(20),3080–3082(2009).4.S.Zhang,D.van der Weide,and J.Olvier,“Superfast phase-shifting method for3-D shape measurement,”Opt.Express18(9),9684–9689(2010).5.G.Sansoni,M.Carocci,and R.Rodella,“Three-dimensional vision based on a combination of gray-code andphase-shift light projection:analysis and compensation of the systematic errors,”Appl.Opt.38(31),6565–6573 (1999).6.S.Zhang,“Flexible3-D shape measurement using projector defocusing:Extended measurement range,”Opt.Lett.35(7),931–933(2010).7.J.Pan,P.S.Huang,and F.Chiang,“Color-coded binary fringe projection technique for3D shape measurement,”Opt.Eng.44(2),023606(2005).8. D.Malacara,ed.,Optical Shop Testing,3rd ed.(John Wiley and Sons,2007).9. D.C.Ghiglia and M.D.Pritt,Two-Dimensional Phase Unwrapping:Theory,Algorithms,and Software(JohnWiley and Sons,1998).10. C.Zhang,P.S.Huang,and F.-P.Chiang,“Microscopic phase-shifting profilometry based on digital micromirrordevice technology,”Appl.Opt.41,5896–5904(2002).1.IntroductionHigh-resolution,real-time3D shape measurement has become a major research subject due to its importance in numerous applications[1].Digital fringe projection techniques have shown their potential for dynamic3D shape measurement.However,when a commercially available #142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8539digital-light-processing(DLP)projector is used,the measurement speed is limited to120Hz, which is its highest pattern switching rate[2].To improve the measurement speed,we have proposed the binary defocusing technique[3]. Specifically,this technique is to generate sinusoidal fringe patterns by properly defocusing the1-bit binary structured paring with a8-bit grayscale sinusoidal fringe pattern, the binary structured one requires smaller data size to represent,thus allows for faster pattern switching rate.For example,the DLP Discovery D4000can switch images at faster than32kHz for1-bit binary patterns,whilst only291Hz for8-bit grayscale ones.If the binary defocusing and a three-step phase-shifting algorithm are used,the achievable3D shape measurement speed is faster than10kHz.In contrast,the same DLP Discovery can only achieve97Hz3D shape measurement speed if8-bit sinusoidal fringe patterns are used.With the binary defocusing technique,we have successfully achieved3D shape measurement at an unprecedentedly high speed:667Hz[4].Instead of10kHz,667Hz speed was limited by the projection source intensity rather than by the capability of the projector or the camera.The superfast3D shape measurement system[4],however,cannot be used to measure mul-tiple moving objects because a single-wavelength three-step phase-shifting algorithm is used. Adding additional coded patterns for fringe order determination has the potential to resolve this problem because it can obtain absolute phase point-by-point[5].But,because the defocusing technique is used,the challenge becomes more complicated.For example,Zhang has demon-strated that by combing with simple binary coding with phase shifting,the point-by-point3D shape measurement could be realized[6].However,because the binary coding technique used was not based on gray code,the decoding process has significant errors due to defocusing.To address this problem,we adopt the gray code[7]in this research.The gray code changes only one bit at one position.However,we found even the gray code has problems for rapid motion capture.This is because when the object is moving rapidly,there might be significant displace-ment between thefirst captured pattern and the last one.Thus the precisely designed code will have misalignment,which may introduce severe errors in thefinal reconstructed results.In this paper,we propose a new phase unwrapping strategy that combines the conventional spatial phase unwrapping algorithm with the gray code to solve the motion induced phase un-wrapping problems.The gray code isfirst used to segment the whole image into separate re-gions.Within each region,the phase differences between neighboring pixels are no more than 2π.Then a conventional unwrapping algorithm is applied to unwrap the phase for each region to obtain relative phase.Finally,the gray code is again used to obtain absolute phase by prop-erly shifting the relative phase for each region by adding multiples of2π.By this means,the phase unwrapping problems caused by motion could be resolved.To verify the performance of the proposed method,a system with a speed of5,000Hz was developed to capture the collision process of two balls.Section2explains the principle of the absolute phase retrieval by combining gray coding with phase shifting.Section3presents the proposed framework to solve the motion induced phase unwrapping problems.Section4shows experimental results,andfinally Section5summarizes this paper.2.Principle of phase-shifting technique2.1.Three-step phase-shifting algorithmPhase-shifting methods are widely used in optical metrology because of their speed and ac-curacy[8].We use a three-step phase-shifting algorithm tofind the phase value.Three fringe #142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8540images with a phase shift of2π/3can be written as:I1(x,y)=I (x,y)+I (x,y)cos(φ−2π/3),(1)I2(x,y)=I (x,y)+I (x,y)cos(φ),(2)I3(x,y)=I (x,y)+I (x,y)cos(φ+2π/3).(3) where I (x,y)is the average intensity,I (x,y)the intensity modulation,andφ(x,y)the phase to be solved for.The phase can be solved for from these equations asφ(x,y)=tan−1 √3(I1−I3)/(2I2−I1−I3).(4)This equation provides the wrapped phase ranging from−πto+πwith2πdiscontinuities. The2πphase jumps can be removed to obtain a continuous phase map by adopting a phase unwrapping algorithm[9].Specifically,the phase unwrapping is to determine fringe order, integer k(x,y),so thatΦ(x,y)=φ(x,y)+k(x,y)×2π.(5) Here,Φ(x,y)is the unwrapped absolute phase.However,such a single three-step phase-shifting algorithm cannot be used to measure multiple objects at a time.2.2.Fringe order k(x,y)determination with gray codeTo simultaneously measure multiple objects,a binary coding method can be used for absolute phase extraction.Integer number k(x,y)in Eq.(5)can be uniquely defined by a sequence of binary(1-bit)patterns.Figure1illustrates two coding examples.It can be seen that for each fringe period,a unique bit sequence(e.g.,0001,0011)is used.However,due to the projector defocusing and discrete camera sampling,the binary state changes in coded patterns become blurred,and this makes the code boundaries difficult to be determined accurately.This problem will be worse if the state changes at the same place for different coded patterns,as illustrated in Fig.1(a):all four patterns change from code0111to1000.To reduce this problem,we adopt a gray code[7]as illustrated in Fig.1(b).Gray code is to ensure that at each point there is just one bit change for all the coded ing this new coding method,the codeword-detecting errors could be reduced and the unwrapping result can be improved.However,when the objects are moving fast,the motion could cause additional phase unwrapping problems even if the gray coding patterns are used.There will be misalign-ment between the extracted wrapped phase and precisely designed codeword.Therefore,severe errors will arise near the jump edges of codeword,which will be illustrated in the Sec.4.(,)x yICode(a)Traditional coding(b)Gray codingFig.1.Gray Vs non-gray coding for point-by-point phase unwrapping.2.3.Phase-to-depth conversionAfter obtaining the absolute phaseΦ(x,y),the depth distribution can be retrieved by adopting a phase-to-height conversion algorithm.In this paper,a reference-plane-based method is adopted#142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8541to convert absolute phase to depth z[10].In brief,assumingΦ(x,y)is the object phase map, andΦrp(x,y)is the reference phase map that is obtained by measuring a uniformflat planar surface,the relationship between the surface height relative to the reference plane and their phase difference,ΔΦ(x,y)=Φ(x,y)−Φrp(x,y),is approximately proportional,that isΔz∝ΔΦ(x,y)=Φ(x,y)−Φrp(x,y).(6)Assuming the reference plane has a depth of z0=0,the absolute depth value for each camera pixel can be represented asz(x,y)=z0+c0×[Φ(x,y)−Φrp(x,y)]=c0×[Φ(x,y)−Φrp(x,y)],(7)where c0is a constant that can be determined through calibration.3.Novel phase unwrapping frameworkTo tackle the phase unwrapping problems caused by motion,a novel phase unwrapping frame-work is proposed.Specifically,we combine the conventional spatial phase unwrapping method with the gray coding method.The framework can be implemented in the following steps: Step1:Use the codeword to segment the whole image into separated regions.For each point, we can detect the codeword jump by comparing it with neighboring points.For each region,the codeword does not change more than1,which means that the phase difference between neigh-bor pixels should not be more than2π.If the jump of one point is more than1,meaning that the phase difference of the point might be more than2π,this point is treated as a discontinuous edge.The entire image is then divided into continuous regions separated by the edges.Step2:Adopt a conventional phase unwrapping algorithm to unwrap the phase for each region to obtain the relative phase map.When the phase map is continuous in one region,a conventional unwrapping algorithm can be applied to obtain the continuous phase map.Be-cause each region is unwrapped individually relative to one point in that region,the unwrapped phase obtained in this step is called relative phase,Φr(x,y).Step3:Use the codeword to adjust the relative phase maps to be absolute ones for all re-gions.Step2only provides relative phase for each region.In order to measure multiple objects simultaneously,absolute phase is needed.To do so,wefirst use the conventional approach to extract the absolute phase mapΦ0(x,y).For each region,the difference between the ab-solute phaseΦ0(x,y)and the relative phaseΦr(x,y)obtained in Step2should be a constant,ΔΦ(x,y)=Φr(x,y)−Φ0(x,y)=Δk(x,y)×2π.Here,Δk(x,y)is an integer and constant for each region.The difference values should be nearly the same for each region except for a few noisy points.To obtainΔk(x,y),we canfirst average all the phase difference values in one region after applying a low-passfilter,and then divide the average value by2πand round it to an integer for each region.OnceΔk(x,y)is obtained,the relative phase can be converted to absolute one by addingΔk(x,y)×2π.Step4:Combine different region phase maps into a single one and convert thefinal un-wrapped phase map to3D shape.This step will generate thefinal complete absolute phase map Φ(x,y)by combing all regions together.Once the absolute phase map is obtained,3D shape can be recovered if the system is calibrated.4.Experimental resultsWe developed a superfast system that is composed of a DLP Discovery projection system,a high-speed CMOS camera,and a self-developed synchronization circuit.The DLP Discovery projection system includes a DLP Discovery board(D4000)(Texas Instruments)with an ALP#142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8542High Speed(Digital Light Innovations).The camera used in this system is Phantom V9.1(Vi-sion Research).The synchronization circuit takes the projection timing signal and sends the trigger signal to the camera for simultaneous image acquisition.Experiments were conducted to capture the collision of two balls.For each3D frame,four binary coded patterns and three phase-shifted fringe patterns are projected.In this research, the patterns are continuously captured at5,000Hz with a480×480image resolution.The exposure time used is200μs.In the experiment,both balls are suspended by strings,and initially in contact with one another.Then,one ball is pulled to a certain height and released to swing like a pendulum until collision occurs with the stationary ball,at the bottom of the pendular arc where the moving ball’s velocity is at its maximum.Figure2(a)shows one fringe pattern and Fig.2(b)shows one of the coded binary patterns.From the four coded binary patterns,the codeword k(x,y),which represents the fringe order of each point,can be extracted as shown in Fig.2(c).From the three phase-shifted fringe patterns,the wrapped phase map can be obtained as shown in Fig.2(d).The codeword can be used to unwrap the phase point-by-point using Eq.(5).Figure2(e)shows the result.It clearly shows some problems on the unwrapped phase(e.g.,stripes and spikes).To highlight these problems,Fig.2(f)shows a view zoomed into a small area.Thisfigure clearly shows that the problematic areas are more than 2pixels wide.This means that a standardfiltering technique cannot be used to eliminate this problem.Since the object is moving rapidly,the coded patterns and the phase-shifted fringe patterns are captured sequentially,and the time lag between thefirst frame and the last frame is1.4ms,the displacement between thefirst captured fringe pattern and the last one could be significant.Thus the precisely designed codeword and the wrapped phase do not match,and the phase cannot be unwrapped correctly by the coded patterns.(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Fig.2.(a)One fringe pattern;(b)One binary pattern;(c)Fringe order k(x,y)extracted fromfour coded patterns;(d)Wrapped phaseφ(x,y)extracted from three phase-shifted fringepatterns;(e)Absolute phaseΦ0(x,y);(f)Zoom-in view of the area within the white windowof(e);(g)Codeword change map;(h)Unwrapped phase for one region by a conventionalunwrapping algorithm;(i)Unwrapped relative phase of all regionsΦr(x,y);(j)Differencebetween relative phase map in(i)and absolute phase map obtained traditionally(e);(k)Absolute phaseΦ(x,y)after Step3;(l)The corresponding area shown in(f).Now we apply the novel phase unwrapping framework we proposed in Sec.3.Figure2(g) shows the map of codeword changes which can be used to segment the images into different regions.By applying a conventional phase unwrapping algorithm to each region,the continuous phase map can be obtained.One of the regions is illustrated in Fig.2(h).However,because#142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8543each region is unwrapped relative to one point within that region,we only obtain relative phase instead of absolute one.For example,Fig.2(i)shows the unwrapped phase for two regions,it clearly shows that they are not absolute because the phase is not monotonic along x-axis.The relative phase mapΦr(x,y)can be converted to absolute one by referring to the absolute phase mapΦ0(x,y)obtained by using a conventional approach as shown in Fig.2(e).Figure2(j)shows the difference phase mapΔΦ(x,y)=Φr(x,y)−Φ0(x,y).Thisfigure indicates that the difference values indeed are nearly the same for each region except for a few noisy points. Once the constants for each region is determined,the relative phase can be converted to be absolute one as show in Fig.2(k).Because the relative phase map for each region is obtained through a conventional phase unwrapping algorithm,it should be smooth.Figure2(l)shows the corresponding area of the absolute phase as shown in Fig.2(f).It can be seen that the noisy areas are gone,and the phase map is smooth.(a)(b)(c)Fig.3.(a)Photograph of one frame;(b)3D result using the conventional method;(c)3Dresult by the proposed phase unwrapping framework.(Media1,Media2,Media3,andMedia4).Figure3shows the experimental result for one frame.Figures3(a),3(b)and3(c)respectively show the photograph of the frame,the3D reconstructed results using the conventional and the proposed technique.It clearly indicates that the spike and large noisy points disappear if the new method is adopted.This means that the new phase unwrapping method can resolve the phase unwrapping problem caused by motion.Several videos are also submitted together with this paper.Media1shows the full speed video(5000/7=714Hz)of the collision process.Media2, Media3,and Media4shows the video at60Hz viewing from different angles.Our experiment shows that the collision process of two balls are well captured and3D information are correctly recovered by adopting the new phase unwrapping strategy and thus verifies the success of the proposed method.5.ConclusionThis paper has presented a novel and efficient phase unwrapping strategy to solve phase un-wrapping problems associated with3D shape measurement of moving objects.Experimental results demonstrate that the proposed method can overcome motion-induced phase unwrapping problems for simultaneously measuring multiple rapidly moving objects.However,the pro-posed method is not trouble free.The assumption of the framework is that within each region there are no step-height changes betweenπand2π.Although this is not likely to happen in most measurement scenarios,if it becomes a problem for some particular applications,various image processing means could be adopted to address it.It is important to note that this proposed method only solves for the phase unwrapping error caused by motion,but cannot reduce the inherent phase error caused by motion.For this system with5,000Hz capture speed and200μs exposure time,we assume the object is motionless#142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8544during three phase-shifted fringe images capture,i.e.,the phase error caused by600μs time lag will not be alleviated.Even with this limitation,the proposed system has achieved a3D shape measurement with an equivalent shutter speed of1,667Hz or600μs shutter time,which can be used to capture relatively rapidly moving objects.#142632 - $15.00 USD Received 17 Feb 2011; revised 30 Mar 2011; accepted 30 Mar 2011; published 18 Apr 2011 (C) 2011 OSA25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8545。

Shape Google:geometric words and expressions for invariant shape retrievalALEXANDER M.BRONSTEINDepartment of Electrical Engineering,T el-Aviv UniversityMICHAEL M.BRONSTEINInstitute of Computational Science,Faculty of Informatics,Universit`a della Svizzera Italiana LEONIDAS J.GUIBASDepartment of Computer Science,Stanford Universityand MAKS OVSJANIKOVInstitute for Computational and Mathematical Engineering,Stanford UniversityThe computer vision and pattern recognition communities have re-cently witnessed a surge of feature-based methods in object recog-nition and image retrieval applications.These methods allow repre-senting images as collections of“visual words”and treat them using text search approaches following the“bag of features”paradigm.In this paper,we explore analogous approaches in the3D world applied to the problem of non-rigid shape retrieval in large ing multiscale diffusion heat kernels as“geometric words”,we construct compact and informative shape descriptors by means of the“bag of features”approach.We also show that considering pairs of“geometric words”(“geometric expressions”)allows creating spatially-sensitive bags of features with better discriminative power.Finally,adopting metric learning approaches,we show that shapes can be efficiently represented as binary codes.Our approach achieves state-of-the-art results on the SHREC2010large-scale shape retrieval benchmark. Categories and Subject Descriptors:1.INTRODUCTIONThe availability of large public-domain databases of3D mod-els such as Google3D Warehouse has created the demand for shape search and retrieval algorithms capable offinding sim-ilar shapes in the same way a search engine responds to text queries.However,while text search methods are sufficiently developed to be ubiquitously used e.g.in a Web application, the search and retrieval of3D shapes remains a challenging problem.Shape retrieval based on text metadata(annotations and tags added by humans)is often not capable of providing the same experience as a text search engine[Min et al.2004]. Content-based shape retrieval using the shape itself as a query and based on the comparison of the geometric and topo-logical properties of shapes is complicated by the fact that many3D objects manifest rich variability,and shape retrieval must often be invariant under different classes of transforma-tions.A particularly challenging setting,which we address in this paper,is the case of non-rigid or deformable shapes, which includes a wide range of shape transformations such as bending and articulated motion.An analogous problem in the image domain is image re-trieval:the problem offinding images depicting similar scenes or objects.Similar to three-dimensional shapes,images may manifest significant variability(Figure1),and the aim of a successful retrieval approach is to be insensitive to such changes,while maintaining high discriminative power.Signif-icant advances have been made in designing efficient image retrieval techniques(see an overview in[Veltkamp and Hage-doorn2001]),but the majority of2D retrieval methods do not immediately generalize to3D shape retrieval[Tangelder and Veltkamp2008].Recently,feature-based methods have gained popularity in the computer vision and pattern recognition communities with the introduction of scale invariant feature transform(SIFT) [Lowe2004]and similar algorithms[Matas et al.2004;Bay et al.2006].The ability of these methods to demonstrate suf-ficiently good performance in many settings,including object recognition and image retrieval and the public availability of the code made SIFT-like approaches a commodity and a de facto standard in a variety of image analysis tasks.One of the strengths of feature-based approaches in image retrieval is that they allow one to think of an image as a col-lection of primitive elements(visual“words”),and use the well-developed methods from text search such as the“bag of features”paradigm.One of the best implementations of these ideas is Video Google,a web application for object search in large collections of images and videos developed in Oxford university by Zisserman and collaborators[Sivic and Zisser-man2003;Chum et al.2007],borrowing its name through an analogy with the famous text search engine.Video Google makes use of feature detectors and descriptors to represent an image as a collection of visual words indexed in a“visual vo-cabulary.”Each image is compactly encoded into a vector of frequencies of occurrences of visual words,a representation referred to as a“bag of features.”Images containing similar ACM Transactions on Graphics,V ol.,No.,Article,Publication date:.2•Bronstein et al.visual information tend to have similar bags of features,and thus comparing bags of features allows retrieving similar im-ages.Zisserman et al.showed that employing weighting schemes for bags of features that take into consideration the average occurrence of visual words in the whole database allows for very accurate retrieval[Sivic and Zisserman2003;Chum et al. 2007].Since the comparison of bags of features usually boils down tofinding weighted correlation between vectors,such a method is suitable for indexing and searching very large (Internet-scale)databases of images.In a follow-up work,Grauman et al.[Jain et al.2008] showed that an optimal weighting of bags of features can be learned from examples of similar and dissimilar im-ages using metric learning approaches[Torralba et al.2008]. Shakhnarovich[2005]proposed the similarity-sensitive hash-ing,which regards metric learning as a boosted classification problem.This method appeared very efficient in learning in-variance to transformations in the context of video retrieval [Bronstein et al.2010c].In[Bronstein et al.2010e],an ex-tension of this approach to the multi-modal setting was pre-sented.In[Strecha et al.2010],similarity-sensitive hashing algorithms were applied to local SIFT descriptors to improve the performance of feature matching in wide-baseline stereo reconstruction problems.Behmo et al.[2008]showed that one of the disadvantages of the bag of features approaches is that they lose informa-tion about the spatial location of features in the image,and proposed the commute graph representation,which partially preserves the spatial information.An extension of this work based on the construction of vocabularies of spatial relations between features was proposed in[Bronstein and Bronstein 2010a].The success of feature-based methods in the computer vi-sion community is the main inspiration for the current paper, where we present a similar paradigm for3D shapes.1.1Related works in the shape analysiscommunityShape retrieval is an established research area with many ap-proaches and methods.For a detailed recent review,we re-fer the reader to[Tangelder and Veltkamp2008].In rigid shape retrieval,global descriptors based on volume and area [Zhang and Chen2001],wavelets[Paquet et al.2000],sta-tistical moments[Kazhdan et al.2003;Novotni and Klein 2003;Tal et al.2001],self-similarity(symmetry)[Kazhdan et al.2004],and distance distributions[Osada et al.2002] were used.Methods reducing the3D shape retrieval to image retrieval use2D views[Funkhouser et al.2003;Chen et al. 2003],slices[Jiantao et al.2004],silhouette and contour de-scriptors[Napol´e on et al.2007].Graph-based methods based on skeletons[Sundar et al.2003]and Reeb graphs[Hilaga et al.2001;Biasotti et al.2003;Tung and Schmitt2005]are capable of dealing with deformations e.g.matching articulated shapes.Lipman and Funkhouser[2009]proposed the M¨o bius voting scheme for sparse shape matching.Isometric shape deformations werefirst explicitly addressed by Elad and Kimmel[2001;2003].The authors used multi-dimensional scaling(MDS)[Borg and Groenen1997;Bron-stein et al.2006]to construct a representation of the intrin-sic geometry of shapes(captured as a matrix of inter-point geodesic distances and referred to as canonical forms)in a low-dimensional Euclidean space.A moment-based shape de-scriptor[Tal et al.2001]was then applied to obtain a shape signature.The method of Elad and Kimmel was extended in [M´e moli and Sapiro2005;Bronstein et al.2006b;2006a]to compare shapes as metric spaces using the Gromov-Hausdorff distance[Gromov1981],which tries tofind the minimum-distortion correspondence between two shapes.Numerically, the Gromov-Hausdorff distance computation can be carried out using a method similar to MDS[Bronstein et al.2006b] or graph labeling[Torresani et al.2008;Wang et al.2010]. In follow-up works,extensions of the Gromov-Hausdorff dis-tance used diffusion geometries[Bronstein et al.2010d],dif-ferent distortion criteria[M´e moli2007;M´e moli2009],local photometric[Thorstensen and Keriven2009]and geometric [Dubrovina and Kimmel2010;Wang et al.2010]data,and third-order[Zeng et al.2010]distortion terms.While allowing for very accurate and theoretically isometry-invariant shape comparison,the main drawback of the Gromov-Hausdorff framework is that it is based on computationally-expensive op-timization.As a result,such methods are mostly suitable for one-to-one or one-to-few shape retrieval cases,even though there have been recent attempts to overcome this difficulty us-ing Gromov-Hausdorff stable invariants[Chazal et al.2009], which can be computed efficiently in practice.Reuter et al.[2005;2009]proposed using Laplace-Beltrami spectra(eigenvalues)as isometry-invariant shape descrip-tors.The authors noted that such descriptors are invariant to isospectral shape transformations,a family theoretically larger than isometries.Rustamov[2007]used an isometry-invariant shape embedding similar to the eigenmaps proposed in[B´e rard et al.1994;Belkin and Niyogi2003;Lafon2004] to create shape representations in the Euclidean space simi-lar in spirit to canonical forms.He used histograms of Eu-clidean distances(using the approach of Osada et al.[Os-ada et al.2002])in the embedding space to compare the shapes.Shape similarity based on the comparison of his-tograms of diffusion distances were used as by Mahmoudi and Sapiro[2009].In[Bronstein et al.2010d;Bronstein and Bron-stein2009;2010b],an intimate relation between this method and the methods of Rustamov[Rustamov2007]was shown, and[Bronstein and Bronstein2010b]showed a more general framework of which[Mahmoudi and Sapiro2009;Rustamov 2007]are particular cases.M´e moli[2009]formulated a proper metric between shapes based on a spectral variant of the Gromov-Wasserstein dis-tance,and showed that shape similarity measures of Reuter et al.[Reuter et al.2005],Rustamov[2007],and Sun et al.ACM Transactions on Graphics,V ol.,No.,Article,Publication date:.Shape Google:geometric words and expressions for invariant shape retrieval•3[2009]can be viewed as a hierarchy of lower bounds of this metric.Local feature-based methods are less common in the shape analysis community than in computer vision,as there is noth-ing equivalent to a robust feature descriptor like SIFT to be universally adopted.We see a few possible reasons.First, one of the important properties of SIFT is its discriminativity combined with robustness to different image transformations. Compared to images,shapes are usually known to be poorer in features,and thus descriptors are less informative.Secondly, unlike images where invariance is usually limited to affine transformations,the degree of invariance required for shapes is usually much larger and includes non-rigid and topology-changing deformations.Feature-based methods have been explored by Pauly et al.[2003]for non-photorealistic rendering.The authors de-tect multiscale robust contour features(somewhat analogous to edges in images).A similar approach was proposed in [Kolomenkin et al.2009].In[Mitra et al.2010],local fea-tures were used to detect shape self-similarity and grid-like structures.Local moments[Clarenz et al.2004]and vol-ume descriptors[Gelfand et al.2005]have been also pro-posed for rigid shape retrieval and correspondence.Shilane and Funkhauser[2006]use a learnable feature detector and de-scriptor maximizing the likelihood of correct retrieval,based on spherical harmonics and surface area relative to the cen-ter of mass.Mitra et al.[2006]proposed a patch-based shape descriptor and showed its application to shape retrieval and comparison.They also showed an efficient hashing scheme of such descriptors.An approach similar to the one presented in this paper was presented in[Li et al.2006]for rigid shapes. Relatively few feature-based methods are invariant to iso-metric deformations by construction.Raviv et al.[Raviv et al. 2007]and Bronstein et al.[Bronstein et al.2009]used his-tograms of local geodesic distances to construct an isometry-invariant local descriptor.In[Bronstein et al.2009],statisti-cal weighting scheme similar to[Sivic and Zisserman2003; Chum et al.2007]was also used to define point-wise sig-nificance measure in the problem of partial shape matching. Ben-Chen et al.[Ben-Chen et al.2008]use conformal factors as isometry-invariant local descriptors.Zaharescu et al.[Za-harescu et al.2009]proposed a SIFT-like descriptor applied to functions defined on manifolds.Sun et al.[Sun et al.2009]in-troduced deformation-invariant descriptors based on diffusion heat kernels(referred to as heat kernel signatures or HKS). More recently,a scale-invariant modification of these descrip-tors(SI-HKS)were constructed in[Bronstein and Kokkinos 2010]using local scale normalization,and a volumetric heat kernel signature(vHKS)was proposed in[Raviv et al.2010]. These former two approaches(HKS and SI-HKS)are adopted here for their discriminativity and efficient computability. 1.2Main contributionIn this paper,we bring the spirit of feature-based computer vision approaches to the problem of non-rigid shapesearch Fig.1.Examples of invariant image(top)and shape(bottom) retrieval.Shown on the left is a query,and on the right a few examples of desired correct matches retrieved from a large database.Transformations shown in image retrieval are view-point variation,different illumination,background variation, occlusion and partially missing data;in shape retrieval,differ-ent non-rigid shape deformations are shown.and retrieval.By analogy to Zisserman’s group works,we call our method Shape Google.The present paper is an extended version of[Ovsjanikov et al.2009],where the approach was first introduced.While working on this paper,we discovered that the use of bags of features for shape representation has been independently developed by Toldo et al.[2009].In this paper,wefirst show a feature detector and descrip-tor based on heat kernels of the Laplace-Beltrami operator,in-spired by Sun et al.[2009].Descriptors are used to construct a vocabulary of geometric words,distributions over which serve as a representation of a shape.This representation is invariant to isometric deformations,robust under a wide class of per-turbations,and allows one to compare shapes undergoing dif-ferent deformations.Secondly,we show that taking into con-sideration the spatial relations between features in an approach similar to commute graphs[Behmo et al.2008]allows improv-ing the retrieval performance.Finally,adopting metric learn-ing techniques widely used in the computer vision community [Jain et al.2008],we show how to represent shapes as compact binary codes that can be efficiently indexed and compared us-ing the Hamming distance.Figure2depicts aflow diagram of the presented approach. The shape is represented as a collection of local feature de-scriptors(either dense or computed at a set of stable points following an optional stage of feature detection).The descrip-tors are then represented by“geometric words”from a“geo-metric vocabulary”using vector quantization,which produces a shape representation as a bag of geometric words or pairs of words(expressions).Finally,similarity sensitive hashing is ap-plied to the bags of features.We emphasize that the presented approach is generic,and different descriptor and detectors can be used depending on the application demands.ACM Transactions on Graphics,V ol.,No.,Article,Publication date:.4•Bronstein etal. (Featuredetection)Featuredescription Vectorquantization Bag ofwords Bag of expressions Fig.2.Flow of the ShapeGoogle algorithm.The rest of this paper is organized as follows.In Section2, we start with a brief overview of feature-based approaches in computer vision,focusing on methods employed in Video Google.In Section3,we formulate a similar approach for shapes.We show how to detect and describe local geomet-ric features.In Section4,we describe the construction of bags of geometric words.In Section5,we explore metric learning techniques for representing shapes as short binary codes using Hamming embedding.Section6shows experimental results, and Section7concludes the paper.2.BACKGROUND:FEA TURE-BASEDMETHODS IN COMPUTER VISIONThe construction of a feature-based representation of an image typically consists of two stages,feature detection and feature description,often combined into a single algorithm.The main goal of a feature detector is tofind stable points or regions in an image that carry significant information on the one hand and can be repeatedly found in transformed versions of the im-age on the other.Since there is no clear definition of what is a feature,different approaches can be employed.For example, in the SIFT method,feature points are located by looking for local maxima of the discrete image Laplacian(approximated as a difference of Gaussians)at different scales.SIFT uses lin-ear scale-space in order to search for feature points that ap-pear at multiple resolutions of the image,which also makes the method scale-invariant[Lowe2004].Maximum stable ex-tremal region(MSER)algorithmfinds level sets in the image which exhibit the smallest variation of area when traversing the level-set graph[Matas et al.2004;Kimmel et al.2010]. Finally,it is possible to select all the points in the image or a regular subsampling thereof as the set of features(in the latter case,the detector is usually referred to as dense[Tola et al. 2008]).The next stage is feature description.A feature descriptor uses a representation of local image information in the neigh-borhood of each feature point.For example,SIFT assigns a 128-dimensional descriptor vector constructed as local his-tograms of image gradient orientations around the point.The descriptor itself is oriented by the dominant gradient direction, which makes it rotation-invariant[Lowe2004].A similar ap-proach,Speeded Up Robust Feature(SURF)transform[Bay et al.2006],uses a64-dimensional descriptor,computedeffi-Fig.3.Representation of text(left)and images(right)using the bags of features paradigm.ciently using integral images.At this stage,the image can be compactly represented by specifying the spatial coordinates of the detected feature points together with the corresponding de-scriptors,which can be presented as vectors.This information allows,for example,finding correspondence between images by matching their descriptors[Lowe2004].In order to reduce the representation size,a vocabulary is constructed by performing vector quantization in the descrip-tor space.Descriptors can be replaced by indices in the vo-cabulary representing visual“words”.Typical vocabulary size can vary from a few thousand[Sivic and Zisserman2003]up to one million words[Chum et al.2007].Aggregating all the indices into a histogram by counting the frequency of appear-ance of each visual word,the bag of features(sometimes also called bag of visual terms or bag of visterms)is constructed (Figure3).After the feature detection and description stages,two im-ages can be compared this way by comparing their bags of features.This way,the image similarity problem is reduced to the problem of comparing vectors of feature frequencies. Typically,weighted correlation or weighted Euclidean dis-tance is used to measure similarity of bags of features.The weights can be chosen in such a way that features frequent in the query shape(high term frequency)and infrequent in the entire database(low document frequency)are assigned aACM Transactions on Graphics,V ol.,No.,Article,Publication date:.Shape Google:geometric words and expressions for invariant shape retrieval•5large weight.The weight is expressed as the ratio of the term frequency and the document frequency(referred to as term-frequency inverse document frequency or tf-idf in search en-gine literature).It was shown in[Sivic and Zisserman2003; Chum et al.2007]that this type of weighted distance is supe-rior to simple,non-weighted,approach.In[Jain et al.2008], it was shown that an optimal weighted distance between bags of features on a given database can be learned by supervised learning from examples of similar and dissimilar images.A few recent papers tried to extend bags of features by tak-ing into consideration spatial information about the features. Marszalek and Schmid[2006]used spatial weighting to re-duce the influence of background clutter(a similar approach was proposed in[Leibe et al.2004]).Grauman and Darrell [2005]proposed comparing distributions of local features us-ing earth mover’s distance(EMD)[Rubner et al.2000],which incorporates spatial distances.In[Lazebnik et al.2006],the spatial structure of features was captured using a multiscale bag of features construction.The representation proposed in [Amores et al.2007]used spatial relations between parts. Behmo et al.[2008]proposed a generalization of bags of features that takes into consideration the spatial relations be-tween the features in the image.In this approach,following the stage of feature detection and description,a feature graph is constructed.The connectivity of the graph is determined by the spatial distance and visual similarity of the features(spatially-and visually-close features are connected).Next,the graph is collapsed by grouping together features whose descriptors are quantized to the same index in the visual vocabulary.The con-nections in the collapsed graph represent the commute time be-tween the graph nodes.This graph can be considered an exten-sion of a bag of features,where there is additional information about the relations between the visual words.In[Bronstein and Bronstein2010a],images were described as histograms of pairs of features and the spatial relations between them(re-ferred to as visual expressions),using a visual vocabulary anda vocabulary of spatial relations.3.LOCAL FEA TURES IN SHAPESTrying to adapt feature-based approached to3D shapes,one needs to have the following in mind.First,the type of invari-ance in non-rigid shapes is different from one required in im-ages.Typically,feature detectors and descriptors in images are made invariant to affine transformations,which account for different possible views of an object captured in the image. In case of non-rigid shapes,the richness of transformations is much larger,including changes in pose,bending,and connec-tivity changes.Since many natural shape deformations such as articulated motion can be approximated by isometries,bas-ing the shape descriptors on intrinsic properties of the shape will make it invariant to such deformations.Second,shapes are typically less rich in features than images,making it harder to detect a large number of stable and repeatable feature points. This poses a challenging tradeoff in feature detection between the number of features required to describe a shape on one hand and the number of features that are repeatable on the other,and motivates our choice to avoid feature detection at all and use dense descriptors instead.Third,unlike images which in the vast majority of applications appear as matrices of pixels,shapes may be often represented as triangular meshes, point clouds,voxels,level sets,etc.Therefore,it is desirable to have local features computable across multiple representa-tions.Finally,since shapes usually do not have a global system of coordinates,the construction of spatial relations between features is a challenging problem.There exists a plethora of local shape detection and descrip-tion algorithms,and below we overview some of them.The reader is referred to the review paper[Bustos et al.2005]and the recent benchmarks[Bronstein et al.2010a;Bronstein et al. 2010b]for additional details.3.1Feature detectorsHarris3D.An effective feature detection method,called the Harris operator,first proposed in images[Harris and Stephens 1988]was extended to3D shapes by Glomb[2009]and Sipi-ran and Bustos[2010].This method is based on measuring variability of the shape in a local neighborhood of the point, byfitting a function to the neighborhood,and identifying fea-ture points as points where the derivatives of this function are high[Bronstein et al.2010a].Mesh DOG.Several methods for feature detection have been inspired by the the difference of Gaussians(DOG),a classical feature detection approach used in computer vision. Zaharescu et al.[2009]introduce the mesh DOG approach byfirst applying Gaussianfiltering to functions(e.g.mean or Gauss curvature)defined on the shape.This creates a repre-sentation of the function in scale space,and feature points are prominent maxima of the scale space across scales.Castellani et al.[2008]apply Gaussianfiltering directly on the mesh ge-ometry,and use a robust method inspired by[Itti et al.1998] to detect feature points as points with greatest displacement in the normal direction.Heat kernel feature detectors.Recently,Sun et al.[2009] and Gebal et al.[2009]introduced feature detection methods based on the heat kernel.These methods define a function on the shape,measuring the amount of heat remaining at a point x after large time t given a point source at x at time0,and de-tect features as local maxima of this function.As these meth-ods are intimately related to our work,we discuss in depth the properties of heat kernels in the following section.3.2Feature descriptorsShape context.Though originally proposed for2D shapes and images[Belongie et al.2002],shape context has also been generalized to3D shapes.For a point x on a shape X,the shape context descriptor is computed as a log-polar histogram of the relative coordinates of the other points(x −x)for all x ∈X.Such a descriptor is translation-invariant and can be made rotation-invariant.It is computable on any kind of shape representation,including point clouds,voxels,and tri-ACM Transactions on Graphics,V ol.,No.,Article,Publication date:.6•Bronstein et al.angular meshes.It is also known to be insensitive to small occlusions and distortion,but in general is not deformation-invariant,which makes it disadvantageous in non-rigid shape analysis applications.Spin images.Perhaps one of the best known classes of fea-ture descriptors are spin images[Johnson and Hebert1999; Andreetto et al.2004;Assfalg et al.2007],which describe the neighborhood of a point byfitting a tangent plane to the surface at the point,and accumulating information about the neighborhood into2D images which can then be directly com-pared.Although these methods can be robust with respect to noise and changes in triangulation,they were originally devel-oped for rigid shape comparison,and are thus very sensitive to non-rigid shape deformations.Mesh HOG.Zaharescu et al.[2009]use the histogram of gradients of a function defined in a local neighborhood of a point,as a point descriptor(similar to the histogram of gradi-ents(HOG)[Dalal and Triggs2005]technique used in com-puter vision).Though Zaharescu et al.show insensitivity of their descriptor to non-rigid deformations,the fact that it is constructed based on k-ring neighborhoods makes it theoreti-cally triangulation-dependent.Heat kernel signatures(HKS).Recently,there has been increased interest in the use of diffusion geometry for shape recognition[Rustamov2007;Ovsjanikov et al.2008;Ma-teus et al.2008;Mahmoudi and Sapiro2009;Bronstein et al. 2010d;Raviv et al.2010].This type of geometry arises from the heat equation,∆X+∂∂tu=0,(1)which governs the conduction of heat u on the surface X(here,∆X denotes the negative semi-definite Laplace-Beltrami op-erator,a generalization of the Laplacian to non-Euclidean do-mains).The fundamental solution K t(x,z)of the heat equa-tion,also called the heat kernel,is the solution of(1)with a point heat source at x(see Figure4).Probabilistically,the heat kernel can also be interpreted as the transition density function of a Brownian motion(continuous analog of a random walk) [Hsu2002;Coifman and Lafon2006;Lafon2004].Sun et al.[Sun et al.2009]proposed using the diagonal of the heat kernel as a local descriptor,referred to as the heat kernel signatures(HKS).For each point x on the shape,its heat kernel signature is an n-dimensional descriptor vector of the formp(x)=c(x)(K t1(x,x),...,K tn(x,x)),(2)where c(x)is chosen in such a way that p(x) 2=1.The HKS descriptor has many advantages,which make it a favorable choice for shape retrieval applications.First,the heat kernel is intrinsic(i.e.,expressible solely in terms of the Riemannian structure of X),and thus invariant under isomet-ric deformations of X.This makes HKS deformation-invariant (Figure5,bottom row left).Second,such a descriptor captures information about the neighborhood of a point x on the shape at a scale defined by t.It captures differential informationin Fig.4.Values of K t(x,x)mapped on the shape(left)and val-ues of K t(x,y)for three different choices of y(marked with black dots in three rightmostfigures).The value t=1024is used.Hotter colors represent smaller values.a small neighborhood of x for small t,and global information about the shape for large values of t.Thus,the n-dimensional feature descriptor vector p(x)can be seen as analogous to the multi-scale feature descriptors used in the computer vi-sion community.Third,for small scales t,the HKS descriptor takes into account local information,which makes topologi-cal noise have only local effect(Figure5,bottom row right). Fourth,Sun et al.prove that if the Laplace-Beltrami operator of a shape is non-degenerate(does not contain repeated eigen-values),then any continuous map that preserves the HKS at every point must be an isometry.This latter property led Sun et al.to call the HKS provably informative.Finally,as will be shown next,the computation of the HKS descriptor relies on the computation of thefirst eigenfunctions and eigenvalues of the Laplace-Beltrami operator,which can be done efficiently and across different shape representations.This makes HKS applicable to different geometric data,though we focus in this paper on shapes represented as triangular meshes.Scale-invariant heat kernel signatures(SI-HKS).A dis-advantage of the HKS is its dependence on the global scale of the shape.If X is globally scaled byβ,the corresponding HKS isβ−2Kβ−2t(x,x).It is possible in theory to perform global normalization of the shape(e.g.normalizing the area or Laplace-Beltrami eigenvalues),but such a normalization is im-possible if the shape has e.g.missing parts.As an alternative, a local normalization was proposed in[Bronstein and Kokki-nos2010]based on the properties of the Fourier transform. By using a logarithmic scale-space t=ατ,global scaling re-sults in HKS amplitude scaling and shift by2logαβin the scale-space.This effect is undone by the following sequence of transformations,p dif(x)=(log Kατ2(x,x)−log Kατ1(x,x),...,log Kατm(x,x)−log Kατm−1(x,x)),ˆp(x)=|(F p dif(x))(ω1,...,ωn)|,(3) where F is the discrete Fourier transform,andω1,...,ωn de-notes a set of frequencies at which the transformed vector is sampled.Taking differences of logarithms removes the scalingACM Transactions on Graphics,V ol.,No.,Article,Publication date:.。