CHINESE JOURNAL OF MECHANICAL ENGINEERING●598●Vol. 25, No. 3, 2012DOI: 10.3901CIME.2012.03.598. available online at www springrink.com; www.cjmenet.com; www.jmenet.com.cnMesh Smoothing Method Based on Local Wave AnalysisQIN Xuja.*, ZHENG Hongbo', CHENG Shiwei', LIU Shishuang', and XU Xiaogang-31 School ofComputer Science and Technology, Zhejiang University of Technology Hangzhou 310023. China2 Department of Equipment system and Automatization, Dalian Naval Academy. Dalian 116018, China3 Slate Key Laboratory of CAD&CG Zhejiang Universit, Hangzhou 310058, ChinaReceived January 27.2011; revised September 6, 2011; acepted Seplember 13, 2011Abstract: As the mesh models usually contain noise data, it is ncessary to eliminate the noises and smooth the mcsh. But existedmethods always lose geometric features during the smoothing process. Hence, the noise is considered as a kind of random signal withhigh frequency, and then the mesh model smoothing is operated with signal processing theory. Local wave analysis is used to deal withgeometric signal, and then a novel mesh smoothing method based on the local wave is proposed. The proposed method includesfllowinlg steps: Fistly, analyze the priciple of local wave decomposition for ID signal. and expand it to 2D signal and 3D sphericalsurface signal processing; Secondly, map the mesh to the spherical surface with parameterization, resample the spherical mesh anddecompose the spherical signals by local wave analysis; Thirdly, propose the coordinate smoothing and radical radius smoothingmethods, the former fiters the mesh points' coordinates by local wave, and the ltter filters the radical radius from their geometric centerto mesh points by local wave; Finally, remove the high-frequency component of spherical signal, and obtain the smooth mesh modelwith inversely mapping from the spherical signal. Several mesh models with Gaussian noise are processed by local wave based methodand other compared methods. The results show that local wave based method can obtain better smoothing performance, and reservemore original geometric features at the same time.Key words: mesh models, local wave analysis, noise, smoothingdiscrete surface's smoothness while maintaining the1 Introductiontopology construction and geometric features of the meshmodels.Triangular mesh is a common representation for discreteAt present, triangular mesh model smoothing methodsmesh models. It is the piecewise linear surface which isare mainly based on the Laplacian smoothing method,ormed by triangles edges and vertexes fromoptimization method"， combination method' ”，andthree-dimension space. Triangular mesh represents thesimple non-iterative method'. Laplacian smoothingsurface of the model with better visual efects, and it alsoethod (51, as the unusually one, defines a Laplaciancontrols the triangulated surfaces amount to meet theoperator to adjust directions for each vertex, and then alongdiferent needs of 3D mesh models. There are importantthe directions move the vertices on a certain speed toapplications in NC machining, rapid prototyping modeling,achieve the purpose of adjusting the mesh. This method canvideo and animation modeling. Triangular mesh is ofteneectively adjust the mesh to the regular shape, getobtained through complex surface of geometric samplingwell-proportioned density and shape. But for the meshpoints by a 3D scanner, and the mesh data of objects ismodel which has uneven distributed mesh and contains agenerated by the topology reconstruction or directlylarge number of irregular triangles, the adjustment ofscanning. Due to various factors, 3D data of mesh modelLaplacian smoothing method will make original modeloften include noises, so as before digital geormetrywith lots of deformations. Then, TAUBIN671 proposed aprocessing, we need mesh smoothing. Smoothing process isweighted Laplacian smoothing method, which controlledto eliminate the local disturbance and the noises of thethe deformations in some degree while efectivelytriangular mesh. The goal of smoothing is to obtain moresuppressed the noises, but generated new disturbance in the●Corresponding author. E-mailqxj@ zjut.edu.cnorigin中國煤化.明proposodThis project is supported by Nanonal Narural Science Foundation ofthe av工yst the mesh verticesChina (Grant No. 61075118. Grant No. 61005056. Grant No. 60975016).accordCNMHGucinthedrctionofNational. Key Tchnlogy Suppont Program of Cbina (Grnt No. vertex hormal vector. Altnougn tnis method obtained better2007BAH11B02). Zhejang_ Provincial Natural Science Foundation ofsmoothing effects, it produced a large number of iregularChna (Grant No. 100880 and Open Projet Program of Staite KeyLaboratory of CAD&CG of China (Grant No. A0906)triangles easily. In order to maintain the border feature,c Cuner Mchanral Engneerng Socnuty and Spnge-Venig Bato Hndelberg 2012CHINESE JOURNAL OF MECHANICAL ENGINEERING●599●FLEISHMAN, et all0I, proposed a bilateral filter method values of envelope formed by the local maximum anfor mesh denoising; JONES, et al". proposed non-iterative minimum are zero .One- dimensional and two dimensionalfeature preserving mesh smoothing method; CHENl12] signal can be decomposed into several IMF by local waveproposed the adaptive smoothing method. In recent years, analysis, and each IMF contains single mode vibrationsthere are other smoothing methods, such as the average only, that means the wave has single components.principal curvatures methodl!), and the grid point distanceequalization methdI4.2.1 Local wave analysis for 1D signalFrom the point of view of signal processing theory, theData is complex and often contains noises in real life.noise is a kind of random high-frequency signal. Its Hence, there is few original data which can be consideredfrequency is higher than an artifially setting threshold, as an IMF directly. We used local wave analysis toand it can be filtered through a variety of spatial and decompose the signal data into several IMF withfrequency domain filters. We can consider the local small multi resolutions. The decomposition process is shown asdisturbance of mesh model as noise, but it has some local follows.coherence in the overall figure. The local small disturbance(1) Initialization:r(1)= x(t),i=1.on the mesh model can be regarded as high frequency(2) Extraction of the IMF i:signal information when compared to the whole mesh1) Initialization h=r(1)， k=l;model. Therefore local wave analysis can be applied to2) Compute the local maximum and minimum valueseliminate high frequency components in the mesh modelof h_(t);(eliminate the noise), and make mesh model with higher3) Connect the maxima and minima of hg ()bylevel smoothness after treatment.spline interpolation, to form the upper envelope andLocal waves are the signal fluctuations at locallower envelope;areas'15-16]For a complex non-stationary signal, it may4) Compute the average of the upper and lowercontain multiple oscillation modes at any time and haveenvelope m. _(1); .5) Compute h(1)= h4-()-m;-();Local wave analysis is a new signal analysis method based6) Compute the termination criterion function:on the empirical mode decomposition!". It decomposesnumbers of instantaneous frequency into different intrinsich .(i)-h()Fmode components at one moment, and equals to the sum ofsD,=(1)these intrinsic mode components signals and a trend item.后h U)From the point view of the basic functions, the basicfunctions of local wave analysis are unfixed andIfSDk >8,(0.2<8< 0.3), thenc(1) = h(1), otherwiseadaptivelb8. The process of local wave analysis can be go to2) and letk=k +1.considered as the adaptive filtering based on the(3) Definition:r(1)= 7-(1)- c(1).narrowband filters.(4) If the extremum points of r() are more than two,This paper frstly introduces the principle of ID signal then go to (2) and let i=i+1; otherwise, end thedecomposition and the local wave analysis of 2D image; decomposition, and here 1(1) is the residual amount of thesecondly transforms spherical signals into image, and original data. Finally, the decomposed results isdescribes the local wave decomposition method for the 3Dspherical signals; thirdly. expounds the local wave analysisx()= 2 Cei()+ rm(1).(2)applications for the mesh model smoothing, and use bothcoordinate smoothing method and radial radius smoothingmethod to smooth the mesh model after the sphericalIn this way, we got n IMF components and a residualparameterization; finally uses both curvature map and component. IMF components represent the original signalsmoothing results map to verify the effects of our methods. that contains characteristics of diferent time scales; theresidual component represents the trend information of the2 Local Wave Analysis for 1D and 2D Signaloriginal signal.Currently, a simulation signal is used to describe theLocal waves are the signal fluctuations at local areas. For local wave analysis and it's shown as follows:a complex non-stationary signal, it may contain multipleoscillation modes at any time and have more than onef()= asin2tf! + bsin2nf:t + csin2nfst.(3)instantaneous frequency at a moment. Local wave analysiscan decompose numbers of instantaneous frequencyl+ Slat a中國煤化Ieludes thee vpes ofmoment into different intrinsic mode functions (IMF). An signal::THC N M H Gsoidal; the other twoIMF has two characters: (1) In the entire data set, the are linear FM signals with varying instantaneous frequency.difference value for the amount of extremum and through0 The parameters a, b, c corresponds to amplitudespoint must be 0 or at most I; (2) For any points, the average respcctively, and f, 2 fs corresponds to the instantaneous●600.QIN Xujia, et al: Mesh Smoothing Method Based on Local Wave Analysisfrequencies, and t represents the time. Ifa=l0 mm, b=5the images, and then constructed the interpolation curvedmm, c=10 mm, and fi=5 Hz, f2=10 Hz,fs=25 Hz, the surface. Denote I(x,y) as the digital images, andsimulation signal and local wave decompositions is shownm(x,y) as the algebra mean of the maximum and minimalin Fig. 1. In order to get clear descriptions, there is onlyenvelope curved surface. There are two methods topart of the data here (sampling frequency is 1 000 Hz). construct envelope curved surface by the local extremumig. 1(a) shows the original signal, Fig. 1(b)-Fig. l(d) points of images: one is based on radial basic functionshow three IMF components of the original signal (RBF) interpolation method, and the other is constructingrespectively, Fig. l(e) shows the trend component of the triangular Bezier curved surface interpolation bylocal wave signal.triangulation from the extreme points. In our research, weused compact supported radial basic function (CSRBF)"9]fast interpolation method to construct extremum envelope2:curved surface.Local wave analysis of 2D is described as follows.i MWMMwWMwMwMSet the difference between image I(x, y) and m(x,y) is0.250.50071.00h(x,y) , and h(x,y) is defined asTime 1/s(a) Original simulation signalh(x,y)= I(x,y)- m(x,y),(4)where h(x,y) is one intermnediate process value of I(x,y).And repeat the above process k times until h(x,y)to be an0250.75intrinsic mode function:Timer/s .(b) IMF1hr(x,y)= h 1)(x,)- mn(x,y).(5)Set F(x,y) as one isolated intrinsic mode function, soi WWWWWWwf(x,y)= he(x,y).(6)(C) IMF2Then separate f(x,y) fom the original data:冒1R(x,y)= (x,y)- f(x,gy).(7)075Denote R(x, y) as the new data, repeat the above processTime1/sn times, and the iterative relationship is shown as follows:(d) IMF3R(x,y)-. f2(x,y)= R(x,y),(8)R (,y)-f(x,y)= R(x,y).The final expression is(e) Trend componentFig. 1. Simulation signal and time domain waveformcomponents decomposition by local wave analysis1(x,)= Sf(x ,)+ R(x,y), .(9)=1As shown in Fig. I, local wave analysis can decomposewhere I(x, y) is the original image data, and f(x, y)complex multi-frequency instantaneous signals into a finiterepresents the small-scale details with the descendingnumber of components, and each component cascales, and R.(x, y) is the final large-scale trend item.accordingly represent the instantaneous frequency well. AsWe took the accumulation-type textures as example ofthe process of decompostion without any filters, functions,local wave analysis, shown in Fig. 2(a) below. The imageand uncertainty principle constraints, it can express the size is 200x200 pixes, and consists of the cosinesystem characteristic information more clearly andcomp中國煤化工wih the sampingaccurately.intervYHCNMH G2.2 Local wave analysis of 2D signal(x,2)= Jur(不,石)+ Jmr(x.*)+ fr(x,x2) (10)Extending ID local wave analysis method to 2D data,such as image, we first selected the local extreme points of whereCHINESE JOURNAL OF MECHANICAL ENGINEERING●601●[r(x, x2) = 40cos(60x) + 40cos(40x2).imprecise IMF at the boundary. As the decomposition{Jmr(x，*)= 50cos(17x)+ 30cos(4x2),(11) levels increase, the boundary errors tend to spread to the[fr(x, 2)= 190cos(0.3x)+ 190cos<(0.2x2).internal areas. There are many related approaches testimate the boundary extremum points [20 21".Fig. 2(a) is the original texture image, and though the 2DIn our surface smoothing treatment, we used two-waylocal wave analysis, IMF 1 and IMF2 are shown in extension approach to unfold spherical signal into planeFig. 2()-Fig. 2(), and Fig. 2(d) is the trend component signal along the spherial. cordinates (p,θ)(∈These three components crrespond to the high frequency [0°, 180 1,0∈[0°, 360 ]). Both the sampling frequency ofintermediate frequency and low frequency respectively for latitude and longitude are 1°, and then the unfolded spherethe original texture image.is one grid in361x181. As shown in Fig. 3, along thedirection of θ to extend 30° on both two sides(θ∈[-30° ,390°])，and the interval [- -30° , 0°] and[330° ,360°] have some overlaps, [360° , 390°] and[0° , 30°]have some overlaps (rctangles flld with pointson the right and left shown in Fig. 3).(-30. -30)(a) Accumulation-typc texture(b) IMFI十-θ(0.0)(180, 300)(c) IMF2(d) Irend component(210. 390)Fig. 2. Decomposition results for the two-dimensionalFig. 3. Plane launched two-way signal extensionimage based on the local wave analysisThen the boundary of grids 421x181 was madeAs shown in Fig. 2, based on the local wave analysis, symmetric extensions: along the both positive and negativehigh frequency components contain the image's details, and directions of φ(φ∈[0°, 30°]U[150° , 180*]) (rectanglesthe trend components contain the main contents of the flled with grids at the top and bottom shown in Fig. 3).image. Hence, if we filter the high frequency components, Finally we got a 2D plane signal in 421x 241 grids. Thewe will smooth the image in some degrees.unfolded plane signal processed by the local wave analysiswill effectively eliminate the boundary effects at the valid2.3 Spherical local wave signal analysissignal parts (ectangles flled with gray at the center shownWe defined the unit sphere by latitude and longitude in Fig 3).grids, which had signal values. Then we sampled the Now, for example, giving the spherical signal withsignals along the spherical coordinate(,日), and set the Gaussian noises as the simulated spherical signals and itsinterval as follows:radius R∈[100, 105.1] pixels, as shown in Fig. 4(a). Weused local wave analysis to process this spherical signal,Oφ=Oθ=1° (φ∈[0°,180°], θ∈[0°,360°]).and the result of smoothing is shown in Fig. 4(b).Then we unfolded the spherical signal to 2D planein 361x181 grids, and used the local wave analysis methodsto process the image. As there was a mapping relationshipbetween the spherical signals and their unfolded planesignals, we could inverse mapping the signals processed bythe local wave analysis back to spherical signals.In the local wave analysis, the efectively processing of中國煤化工the boundary is a key to ensure the accuracy of the signaldecomposition and reconstruction. As the maximum and(a.TYHCNMHG)Smoohing resulminimum of boundary points are estimate values, the mean( iaussian noisescurve at the boundary is imprecise, and it also leads to theFig. 4. Smoothing of spherical signals●602●QIN Xujia, et al: Mcsh Smoothing Method Based on Local Wave AnalysisIn fact, the signal smoothing result which removed the vertex split operations, we use local parametric informationGaussian noises was fitered the highfrequency to place two splited vertexes in the unit sphere. As shown incomponents, while it maintained mesh model topology Fig. 5, when all the vertex split operations are finished,information and geometry characters.spherical parameterization of original mesh M is done.Parametric mesh is Mo. The base mash Mo is a convex3 Spherical Parameterization of Mesh Model polyhedron, and sphere mesh M' is obtained through thecenter projection of Mo.Transform analysis is based on specific parameterdomains. For example, the Fourier transform and wavelet 4 Mesh Smoothing Based on Local Wavetransform of 2D image are all based on plane rectangularAnalysisdomain. Previous geometric signal processing methodsmostly parameterize 3D signal to plane rectangular domain, 4.1 Re-sampling and inverse mapping of sphere meshand then they used 2D transform for further process. As the sphere mesh obtained by spherical parameteri-However, the plane parameterization of 3D signals will zation was iregular, we used the local wave analysis ofproduce more intense deformations and distortions. Hence, spherical signal method (see section 2.3) to deal with thewe used spherical parameterization method to map mesh sphere mesh. The iregular spherical surface should bemodel to the unit sphere, and then processed spherical regularly re-sampled by longitudes and latitudes, and thesignals with local wave analysis.interpolation was used to calculate the geometricOur spherical parameterization algorithm was based on coordinates of re-sampling mesh point corresponding to thethe mesh simplification and spherical recovery algorithm original model.proposed by ZHOU, et al22). The algorithm was divided Through spherical parameterizations, each pointinto two parts: simplification and vertex placement. Q(x,y,z) of original mesh corresponds to the sphere meshAccording to the local parameterizations information of the points Q'(o, 0), where x, y, z are the geometric coordinatessimplification, it added new vertexes on the sphere. Fig.5 and p,θ are the spherical coordinates. Set regularillustrates the two main steps of the mesh spherical re-sampling mesh point P is inside the triangle ABC on theparameterization algorithm.spherical mesh, as shown in Fig. 6.M= M,MEdgecllapsecollapseleVenexVertexsplitcfig.6. Coordinates calculationof the triangular mesh pointInterpolation of barycentric coordinates is used to obtainFig. 5. Progressive spherical parameterization algorithmspherical coordinates and geometric coordinatescorresponding to the original mesh of point P:(1) Generate the progressive mesh representation withlocal parameterization information. Loop in the process ofPB(12)edge collapse simplifcation for the original mesh M(M,)until the simplified mesh into a convex polyhedron (theconvex polyhedron is always existed, and the worst is a where , m, n are the area coordinates of point P inside thetetrahedron), which is called base mesh. For each edge triangle, and the area coordinates are calculated as follows:collapse simplification operation, two vertexes are deleted,and they are parameterized to the simplified meshS(OPAB)respectively. The local parameterization information is中國煤化工SiABO，(3)recorded in the corresponding vertex split operation forfYHCNMHGprogressive mesh representation.where S is the triangle's area.(2) Process the records in vertex split operations iThen we used the coordinate x, y,z of MABC in thereverse order based on the initial spbere mesh. For each original mesh and the area coordinate of point P toCHINESE JOURNAL OF MECHANICAL ENGINEERING●603●calculated the coordinate x,y.z of point P:(3) Decompose spherical signals and fitler out highfrequency components of IMF by local wave analysis ofspherical signal. x',y',z' are regarded as the spherical(14) signals, and then we obtain the final gcometric coordinates=l|ya +myo+nyc(x",y",z")of the re-sampling mesh.(4) Obtain new geometric cordinates (x"y",z") ofSince the mapping relationships are existed between thethe re-sampling mesh points by bilinear interpolation, andthese coordinates are also the coordinates of point P.re-sampling mesh and the parametric spherical mesh,Coordinate smoothing method can efectively eliminateparametric spherical mesh and original mesh, we can easilythe high frequency components of the mesh model,obain the coordinate space of the original model by localimplement model smoothing and save the geometrywave based mesh re-sampling. Bilinear interpolation offeatures of original mesh. Fig. 8 is an example of thelatitude and longitude grid is adopted to get parametric gridcoordinate smoothing method's effect.based on re-sampling mesh. Fig. 7 is an example ofmapping and inverse mapping for mesh mdel.(a)Original model(b) Model afier smoothingFig. 8. Example of coordinate smoothing method(a) Original model(b) Spherical paramcterization4.3 Radial radius smoothing methodThe noise of mesh model is local disturbance, and it maycause local diturbance from the mesh points to thegeometric center. Radial radius smoothing method is tocalculate the distances from the geometric center to eachother points, and then local wave analysis is used to processthe distance value and then achieve the smoothing model.The descriptions of radial radius algorithm are shown as(C) Regular spherical sampling(d) Inverse mapping basedfollows.on re-sampling(1) Calculate the geomeric centerO(xo, yo,zo)of theFig 7. Maping and inverse mapping of sperical mesh modeloriginal model. First we calculate the radius 1 of eachmesh point relative to the center o, and then calculate unitradial vector n,(naeony,ng). The formulas are shown as4.2 Coordinate smoothing methodWe processed each three-dimensional coordinate xyzfollows:of mesh model by local wave analysis method, andoreconected the new obtained thee-dimensional(15)coordinates by triangular relationship of original model.Then we got a new mesh model, and obtained tho)smoothing efct. The descrptions of coordinate smoothingalgorithm are shown as fllows.r= J(x-x0)+(),-o)"+(z-zoY，(6)(1) Map each point P(x,y,z) of mesh model. topoint Q(0,B)in the unit sphere afer. the spherical|x-xo|parameterization of the original mesh, and the geometricny ==y- Y(17)coordinates are still recorded as(x,y,z)●(2) Generate the grids with regular latitudes andlongitudes through the re-sampling of spherical中國煤化工parameterization mesh. Then by the coordinateEach puntra,ye, UI uic HIESH model is mapped to.CHCNMH G of the original mesh.interpolation of center of gravity, we can obtain thegeometric coordinate (x',y',z')of point M(φ',的in the point Q(p,θ)in the unit sphere-point, and the radiallatitude and longitude grids.position vector r' of point P is recorded on theCHINESE JOURNAL OF MECHANICAL ENGINEERING●607●Tel: +86-571-85290385; E- mail: qxj@zjut.edu.cnE-mail: swc@zjut.edu.cnZHENG Hongbo, bom in 1977, is currently a lecturer at Zhejiang LIU Shishuang, bom in 1981, is currently a master candidate atUniversiry of Technology China, and working toward her PhD Zhejiang University of Technology China.degree at East China Normal University, China. She received her E-mail: lavender5308@hotmail.commaster degree from Guanxi Teachers Education University, China,in 2003.XU Xiaogang, bom in 1967, is currently a professor and a PhDE-mail: zhb@zjut.edu.cncandidate supervisor at Dalian Naval Academy, China. HCHENG Shiwei, bor in 1981, is currently a lecturer at Zhejiangreceived his PhD degree from Dalian University of Technologv,University of Technology, China. He received his PhD degreeChina, in 1999. His research interests include computer graphicsfrom Zhejiang University, China, in 2009. His rescarch interestsand signal processing.E-mail: xxgang@cad.zju.edu.cninclude computer graphics and human-computer interaction.中國煤化工MHCNMHG