Obtaining malignant melanoma indicators through statistical analysis of 3D skin surface disruptions

Background/purpose: It has been observed that disruptions in skin patterns are larger for malignant melanoma (MM) than benign lesions. In order to extend the classification results achieved for 2D skin patterns, this work intends to investigate the feasibility of lesion classification using 3D skin surface texture, in the form of surface normals acquired from a previously built six‐light photometric stereo device.

M ALIGNANT MELANOMA (MM) is a life-threatening dermatological disease, the successful treatment of which relies heavily on early detection. Employment of a non-invasive automatic analysis method is a useful approach for early differentiation between melanoma and benign pigmented lesions. A group of 2D skin lesion patterns called ABCD (Asymmetry, Border, Colour and Diameter) rules have been used conventionally for the differentiation between benign lesions and MMs; however, only moderate diagnosis accuracy has been achieved (1). To obtain better diagnosis results, new parameters, which can provide additional information, to indicate the presence of melanoma are needed.
Generally, the topography of the skin is directly related to the cell growth patterns under the skin surface. A periodic updating of the skin cells generates a net-like surface with regular geometrical microstructures appearing as a network of triangles and polygons. The important pigmentation cells, melanocytes, may grow up abnormally when the skin receives excessive ultraviolet light radiation or endures some other extreme environment change. Various important skin conditions may present indicators in the form of significant changes in 2D and 3D texture. As a typical example, MM may pathologically occur as a lesion that is variegated in colour and with distortion of the line or net pattern on the skin. Specifically the variable growth rate of the atypical melanocytes not only causes the melanin to be distributed unevenly, resulting in multiple shades of brown or red-brown colour, but it also makes the in situ lesion exhibit a markedly asymmetrical shape with indentations and protrusions (2). Therefore, in addition to the pigmentation information, which is useful for differentiating MMs from benign moles, geometrical microstructure of the 3D skin surfaces may provide another important clue for evaluation of the skin conditions (3).
By extracting a skin texture feature called skin line pattern from macroscopic skin images, MM and benign moles can be distinguished through  (4,5). Some new features associated with the skin line pattern like the skin line direction (SLD) and the skin line intensity have also been found to be useful and complementary to the ABCD feathers (6). However, the skin line patterns extracted from 2D images can only be considered as a kind of pseudo profile of skin surface, as a single 2D image is not enough to describe a 3D object completely. Additionally, the white light clinical images used in the above work are illuminated by side light(s) and are subject to specularities and shadows. Although employment of multiple illumination processing can reduce the effects, it could also reduce the fine details of skin texture patterns due to an averaging procedure (7).
Acquiring 3D data and performing non-invasive 3D analysis has long been an important subject in the development of medical imaging and analysis techniques. Nevertheless, its involvement in dermatological diagnosis is still at an early stage with the majority of previous research are based on analysis of 2D skin data. This is not only because 2D images still comprise the majority of accessible dermatological information, but because it is difficult to capture 3D topography data of the skin due to the dynamic, variegated and multilayered characteristics of the skin.
Handels et al. (8) collected skin 3D profile data using a high-resolution laser profilometer to scan silicon skin replicas. The traditional 2D image features like co-occurrence matrices, Fourier features and fractal features were used on the 3D data and fused with artificial intelligence algorithms to provide a significant improvement in classification performance. Unfortunately, the long scanning time and high price limit the utility of the laser scanning method. Also their research findings are only based on skin replica rather than in vivo skin.
The authors' previous work (9) has proposed a six-light source photometric stereo device that provides a detailed and comprehensive description of the skin. It can obtain both accurate 2D surface reflectance and surface normal map independent from external environment effects such as ambient illumination. Preliminary evidence suggests that the information obtained should be more suitable than conventional photographs for the purposes of visual and automated assessment of skin.
In contrast to the 2D skin line pattern approach, this work seeks to investigate whether analysis of skin surface normal data can derive discriminating MM indicators, in the form of 3D skin surface disruptions, for the differentiation between MMs and benign lesions. The proposed approach will analyse the skin surface normal data through their tilt and slant directions, the so-called skin tilt pattern and skin slant pattern. The basic idea is to simulate 'pseudo-normal' tilt and slant patterns using a 2D Gaussian function, and compare them with a lesion's observed tilt and slant patterns. The differences associated with these two patterns are estimated as the disruptions in skin tilt pattern and skin slant pattern, respectively. It is expected that for a MM, the estimated disruptions will be larger than for a benign lesion.

Material and Methods
The concept of photometric stereo and some notations The surface normal vectors used for the analysis in the current study are obtained from a six-light photometric stereo system; its theoretical foundations are briefly outlined here. For an ideal Lambertian surface, the image irradiance equation can be expressed as where a and b are tilt and slant angles of the illuminants, p 5 dz/dx and q 5 dz/dy (where z is the surface height) are the x-axis component and y-axis component of the surface gradients, respectively, and r is the surface reflection rate (albedo). Although the surface height can be reconstructed through the combination of p and q, this process is subject to a number of limitations and can introduce errors, so the effort of the proposed approach will be focused on analysing the 3D surface normals directly. From Eq. (1), we have where T denotes the transpose. According to Horn (10) ¼ ðn x ; n y ; n z Þ ¼ ðÀp; Àq; 1Þ where N is the unit surface normal vector n x , n y and n z its components in x-axis, y-axis and z-axis, respectively. Hence, the surface gradient vectors can also be calculated as Because a unit surface normal has the modular of 1, only three variables are unknown including the albedo and any two from n x , n y and n z . So we need at least three lights to solve the equation. All other lights offer redundant information which can be used for detecting problematic pixels under specular and shadows and allow them to be removed from the computation of the skin surface normal and reflectance images (11). Figure 1 shows that every surface normal on the surface can be identified by its tilt direction f and slant direction y, which can also be represented by surface gradients as where G is the projected vector of the surface normal vector N on the image plane. It can be observed that Eq. (6) is made up of the surface gradient direction and Eq. (7) is related to the surface gradient magnitude. The employment of the tilt direction f and slant direction y will give more instinctive physical meanings than that of surface gradient or normal. Meanwhile, analysing those two quantities is the equivalent of analysing the surface normals themselves directly. Therefore, the surface normals in the tilt and slant direction may be referred as skin tilt pattern and skin slant pattern.

3D skin texture
It should be noted here that the concept of surface normals derived from the six-light photometric stereo should not be confused with 2D image gradients. This is because surface normals belong to surface topographic texture (or 3D texture) (12). It is quite different from 2D image texture including the image gradients, which inherently are the features formed by pixel intensities under a particular illumination direction, and so is illumination dependent. This is also the main difference between the 2D skin line pattern and the proposed 3D skin texture. The examples for a benign seborrhoeic keratosis and a MM are shown in Fig. 2 where 2D skin line pattern are represented as flows of SLDs, and 3D skin texture profiles are represented as needle maps of surface gradients, which have a trigonometric relation with surface normals as described in Eq. (4).

A variable 3D skin model
Having obtained the 3D surface normal vectors for each lesion, the next step is to estimate their respective skin disruptions. In order to achieve this, a reference skin model needs to be established, whose surface normals can be used as a reference to be compared with the observed surface normals of the lesion.
Depending on lesion class, a lesion's topography can differ from near-flat or slightly raised (such as a benign compound nevus or a junctional nevus) to non-flat or elevated (such as a dome-shaped intradermal nevus) shape. For lesions with flat topography, a reference skin model may be established as a map of constant surface normals. The amount of surface disruptions can be measured directly by calculating surface normal differences between a lesion and the reference skin in the tilt and slant pattern. This approach is feasible for measuring flat lesions but not for raised lesions with a cross-section tended to a hemisphere. For raised lesions, the underlying surface topography will influence the surface normal distributions. As a result of this, previously used surface normal patterns for textured skin are no longer valid. Figure 3 illustrates this problem using a texture composed of regular triangles where the textured skin patterns appear as isosceles, with a disrupted skin pattern shown as a dotted scalene triangle on the left of the figure. On the right of the figure, the dotted hemisphere represents a hemispherical-shaped lesion, on which the skin texture patterns are regularly distanced and well maintained. Here, the surface normal vector N bs from the hemispherical-shaped lesion has the same direction as the one N ds from the disrupted flat skin pattern. If both are compared with the surface normal N s from the reference skin model, they will be recognised as the same surface normal and have the same degree of directional disruption. However, this clearly is not the case as one is from an area of simulated benign skin and the other is from an area of simulated malignant skin.
So in this case, a scheme which represents the topography of the reference skin model that variable and which can be chosen adaptively according to each lesion's topographical characteristics seems reasonable. A natural candidate in this category would be a 2D Gaussian function. The reasons for this are, firstly, it is the most frequent distribution in real life and is widely  used in various parametric statistical hypothesis and analyses; secondly, its flexibility and variability allows it to approximate a wide range of 3D topographies including topographies with sharp protrusions where the variances are set small and the amplitude large, near-flat topographies where the variances are large and the amplitude low, and hemispherical topographies where only the central part of the Gaussian envelope is used. Thirdly, the symmetrical shape of the Gaussian envelope allows asymmetry analysis of surface normals. Hence, now the question in finding this reference model becomes optimising the parameters of the Gaussian that best fits a lesion's observed surface normals.

Disruptions in skin tilt pattern
The shape of a Gaussian distribution may be denoted as where A* is the amplitude, ðu Ã c ; v Ã c Þ is the centre and s* is the variance of the Gaussian function (which is equal in both x-axis and y-axis).
As shown in Eq. (6), the skin tilt pattern is a direction defined entirely within the lesion area on the image plane, so it can be used to estimate the centre location of a Gaussian distribution. By projecting the Gaussian envelope onto the image plane, it is an isotropic distribution of tilt patterns centred at ðu Ã c ; v Ã c Þ. Regarding an isotropic distribution, as shown by Fig. 4, the tilt pattern for every point (u,v) of the isotropic distribution can be determined by its distance to the distribution centre ðu Ã c ; v Ã c Þ in the x-axis and y-axis. Hence, by specifying the centre of the Gaussian function, a map of the tilt patterns for every point in the image can be generated. Because the centre of Gaussian function can be anywhere on the image, we can generate a series of Gaussian envelopes whose number depends on the number of points within the image. In relation to the skin reference model, the objective is to find a surface description that is closest to a lesion, so the centre of the Gaussian function will have to be the one whose Gaussian function best fits (or more formally, giving the minimum directional differences) a lesion's acquired surface normal distribution from the Skin Analyser. This can be described by the following equation where ðu Ã c ; v Ã c Þ is the estimated centre of the Gaussian function where the star sign '*' denotes an estimated variable, S l denotes the lesion region, ||.|| denotes the Euclidean distance. f is the acquired tilt direction and fÃis the tilt pattern of the Gaussian function. Upon finding this distribution centre, the associated directional differences are defined as the Overall disruptions in skin Tilt pattern (OT) between the skin tilt patterns f Ã min of the best fit Gaussian function and the existing skin tilt patterns f Disruptions in skin slant pattern So far, for finding the centre location of the Gaussian distribution, all the computations are limited to the x-y plane and the tilt pattern. To determine the exact topography of the Gaussian distribution, the slant pattern should also be used.
Because the topography of a Gaussian distribution is dependant on it variance and amplitude, by varying the variance and the amplitude, different Gaussian distributions with different topographies are generated. Among them, only the topography or portions of the topography whose skin slant patterns best fit a lesion's acquired skin slant patterns are used to quantify the surface disruptions. This can be defined by the following equation: ðA Ã ; s Ã Þ ¼ arg min X ðu;vÞ2S l y Ã ðu; vÞ À yðu; vÞ k kð11Þ where y ¼ cos À1 n z ð Þ ¼ cos À1 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 þ q 2 þ 1 p ! Fig. 4. The tilt direction at a point (u,v) can be determined by its distance to the distribution centre ðu Ã c ; v Ã c Þ in x-axis and y-axis, hence it can be calculated as fÃ ¼ tan À1 ðjjv À v Ã c jj=jju À u Ã c jjÞ where ||.|| denotes the Euclidean distance. and where p and q are the acquired surface gradients (trigonometric variants of 3D surface normals) in the x-axis and y-axis of a lesion from the skin analyser and (Z Ã x ; Z Ã y ) are the surface gradients in the x-axis and y-axis of the Gaussian functions. The Levenberg-Marquardt method (13) was used to estimate the best fit parameters that solve the non-linear optimisation problem in Eq. (11).
Upon determination of s* and the resultant Gaussian topography or portions of the Gaussian topography, the associated differences are used to define the degree of disruptions in skin slant pattern. The Overall disruptions in skin Slant pattern (OS) is defined as the sum of the directional differences between the skin slant patterns y Ã min of best fit Gaussian function and the acquired skin slant patterns y: In the following experiment, the proposed two features including the skin tilt pattern feature (or OT in Eq. (10)) and the skin slant pattern feature (or OS in Eq. (12)) will be used for lesion classification. To tackle the false 'disruptions' caused by noise, the skin tilt pattern and skin slant pattern disruptions are both enhanced by three iterations of median filters with the window size of 3 Â 3. The discriminating power of the proposed two pattern features will be compared with the best 2D skin line pattern feature (5,6), SLD over the same sample set. The reason to compare with the 2D SLD is because it is the best single 2D feature (among others which include the classic ABCD features and skin line intensity feature) according to (6).

Results
The experiments are carried out on 39 lesions of which 11 are MMs. The 28 benign lesions consist of eight classes including two dermatofibromas (DF), four intradermal nevi (IN), two hyperkeratotic squamous papillomas (HSP), eight compound nevi (CPN), seven seborrhoeic keratoses (SK)s, two dysplastic nevi (DN), two congenital nevi (CGN) and one junctional nevus (JN). All of these lesions were captured at the collaborating institutes using the skin analyser.
The 28 benign lesions are subsequently split into four sample groups as listed in Table 1, with sample group 1 includes one JN and eight CN based on the fact that both lesions are typically small, smooth and slightly raised. Sample group 2 includes seven SKs which are among the most common classes of benign lesions and have a distinct appearance from other benign lesions, sample group 3 includes four IN and two DFs based on the fact that they both have raised and nodular shape. Sample group 4 is made up of the rest of the benign lesions including two DN, two CGN, two HSPs and one blue nevus.
The scatter plots for the skin tilt pattern and skin slant pattern are shown in Figs 5 and 6, respectively. Using the same data samples, the scatter plot for the SLD based on (5) is illustrated in Fig. 7. The discriminating power of each feature is judged by the area under the receiver operating characteristic (ROC) curve or ROC area. The ROC curves for different groups of benign lesions vs. the same 11 MMs are shown in Fig. 8. The corresponding ROC results are shown in Table 2.
As can be seen from Table 2, the overall ROC areas for skin tilt pattern and skin slant pattern for all classes of lesions are 0.73 and 0.85, respectively, both are better than 0.65 achieved by the SLD. Although this ROC result for SLD is not as high as the one reported by She et al. (6), it should be noted that only two types of benign lesions, i.e. JN and CPN, are used in their experiment. In the current experiment, for sample group 1 which consists of JN and CPN, the result for SLD improves to 0.80; a figure in parallel with the one reported by (6) and is equally best with the skin slant feature. Apart from this sample group, the skin slant pattern feature also shows the best results for sample group 2, sample group 3 and all lesion samples. This is probably because skin slant pattern is defined in the z-axis direction, which can reveal extra dimensional information.
In comparison, SLD is inherently a 2D feature while skin tilt pattern is entirely defined on the 2D x-y plane. It can be seen that from Table 2, the classification results for both the skin slant pattern feature and the SLD for sample group 4 which consists of two DN, two CGN and two HSPs are worse than all the other sample groups, this is not surprising as the high difficulty of classifying DN as benign lesions have been generally acknowledged. In fact, the DN and CGN used in the experiments have even been mistaken by the experienced dermatologists at the collaborating clinics and all have subsequently been sent for biopsies. Still substantial improvements over the skin line pattern technique are achieved by the proposed two pattern features with a promising result of 0.83 for the skin tilt pattern feature and 0.77 for the skin slant pattern feature.

Discussion
Prior work using the skin line pattern technique for the diagnosis of MM has used a conventional 2D approach. The classification performance thereby obtained for separating junctional and CPN from MM have in fact been promising. However, the approach has three weaknesses, which should be considered. Firstly, the classification performance relies heavily on obtaining accurate lesion intensities. According to the nature of lesion reflectance, lighting distortions such as specular reflections and shadows are not uncommon. In many cases, these lighting distortions can affect skin line recovery and the classification performance. Secondly, because the computation method (5) for calculating SLD used signed average, it is subject to cancellingout effects, which will make the classification hard to comprehend. This is evident in the more recent skin line pattern publication (14), where only part of the lesion is used for classification, which potentially can reduce these cancelling-out effects. Thirdly, from the generality point of view, the recovery of skin line patterns is sensitive to the value for the 9 Â 9 high-pass filtering scheme used at the beginning for smoothing the skin photographs. With different value (or bandwidth more formally) for smoothing, it is likely that distinct skin line patterns will be generated that could affect the subsequent classification performance.
The proposed statistical analysis on the 3D skin texture does not have the same problems for the following reasons. Firstly, it takes advantages of the additional light sources and images in the photometric stereo technique to isolate lighting/ viewing distortions and improve classification performance. Secondly, the idea of the proposed approach is to use Gaussian functions to simulate    a series of reference skin models, and to find the reference skin model that best fits (or closest to) the surface normal distribution of a lesion. Euclidean distance is used to measure these directional differences and there is not a cancelling-out problem. Third, parameters used in the proposed 3D texture analysis are not set before doing the experiments, those parameters include the centre ðu Ã c ; v Ã c Þ, the amplitude A and the variance s for Gaussian functions. They are found through two separate optimisation procedures (in skin tilt and slant pattern respectively), which minimise the anisotropic differences between the Gaussian functions and the actual surface normal distributions of a lesion.
However, the fundamental difference between the skin line pattern approach and the proposed approach still lies on the difference between 2D skin texture and 3D skin texture. This is because the former one reflects pixel intensity features, while the latter one reflects the inherent objectivity of a lesion. Therefore the proposed analysis of 3D skin texture is independent of the existing 2D techniques, so it provides a new and distinct method for MM diagnosis.

Conclusions
This paper has demonstrated an important application of 3D skin texture for computer-assisted diagnosis of MM in vivo. By taking advantage of the extra dimensional information, preliminary studies suggest that some improvements over the existing 2D skin line pattern approach for separating MM from benign lesions. In terms of methodology, the proposed approach employs a parametric technique, i.e. Gaussian function to adaptively model the surface topographical characteristics of both raised lesions and flat lesions. Whether a reference skin model derived from non-parametric statistical models can improve the classification performance will be investigated in future work.