The surface roughness of several stylolites in limestones was measured using high resolution laser profilometry. The 1D signals obtained were statistically analyzed to determine the scaling behavior and calculate a roughness exponent, also called Hurst exponent. Statistical methods based on the characterization of a single Hurst exponent imply strong assumptions on the mathematical characteristics of the signal: the derivative of the signal (or local increments) should be stationary and have finite variance. The analysis of the measured stylolites show that these properties are not always verified simultaneously. The stylolite profiles show persistence and jumps and several stylolites are not regular, with alternating regular and irregular portions. A new statistical method is proposed here, based on a non-stationary but Gaussian model, to estimate the roughness of the profiles and quantify the heterogeneity of stylolites. This statistical method is based on two parameters: the local roughness (H) which describes the local amplitude of the stylolite, and the amount of irregularities on the signal (mu), which can be linked to the heterogeneities initially present in the rock before the stylolite formed. Using this technique, a classification of the stylolites in two families is proposed: those for which the morphology is homogeneous everywhere and those with alternating regular and irregular portions.