Application of Fractal and Gray-Level Co-Occurrence Matrix Indices to Assess the Forest Dynamics in the Curvature Carpathians-Romania

: The mountain ecosystems face signiﬁcant damage from deforestation and environmental forest changes. We investigated the evolution of tree types of cover areas, deforested areas and total deforested areas from Curvature Carpathians using Gray-Level Co-occurrence Matrix and fractal analysis. The forest dynamics mapping was one of the main objectives of this study and it was carried out using multiple fractal and GLCM indices. We approached the analysis of satellite forest images by calculation of four fractal indices such as Pyramid dimension, Cube Counting Dimension, Fractal Fragmentation-Compaction Index and Tug-of-War lacunarity. We also calculated fractal dimension because it is an index of complexity comparing how the detail in a pattern changes with the scale at which it is measured. Fractal dimension is useful for estimation of irregularity or roughness of fractal and natural objects that do not conform to Euclidian geometry. While the fractal dimension quantiﬁes how much space is occupied, the Tug-of-War lacunarity complements fractal dimension with its ability to quantify how space is occupied. Analysis was further supplemented by the Gray-Level Co-occurrence Matrix analysis because it quantiﬁes spatial probability distributions of gray level values between pixel pairs within an image. The calculated Gray-Level Co-occurrence Matrix features included Angular Second Moment, Contrast, Correlation, Inverse Di ﬀ erence Moment and Entropy. Such comprehensive analysis has the advantage of combining fractal analysis that extracts quantitative information about the morphological complexity of the image with the spatial distribution of the gray pixel intensities as calculated by the co-occurrence features provided by Gray-Level Co-occurrence Matrix. Evolution of deforested areas, expansion of agricultural land and the increased


Introduction
In Earth ecosystems, human activity influences the biodiversity, climate change and a natural order of ecosystems by inducing new pressure elements such as habitats restraining, increasing the species competitivity and migration because of lack of food resources-nourishment. The frequency of extreme events occurs, once they were considered "exceptional" and relegated to long turnaround times, and they led to reconsider the carrying capacity of rivers and the downsizing of the banks enough to contain the floods no longer [1,2]. The forest resources are reduced due to timber production [3] and legal or illegal practices [1] together with climate change in Europe [4][5][6]. Forest areas have increased over the centuries, but after World War II, in socialist countries from South-Eastern Europe the forest areas dramatically decreased because of a high percentage of state-owned forests in Romania (31.2%) [4].
The political changes in the post-socialist period resulted in rapid modification of regulations [7], which liberalized deforestation. Forest management and large forest harvest in Curvature Carpathians are reflected in the increase of forest cover from 1924 to 2014 from 20.1% to 70% [1]. The annual forest disturbance was lower from 2000-2013 (0-0.5%) than 1912-1922 (0.6-1.5%). Percentage of forest in Romania was 25.47% between 1900-1920 with a predominance of deciduous species (52%), coniferous species (25%) and oaks (23%) [1,8]. The Food and Agriculture Organization of the United Nations (FAO) Report from 1986 confirmed that 66% of the forests in Romania are found in the mountains (30% of the country) [9]. From the 1.529 million ha of non-state forests that represent 23% of the total forest area of Romania in 2010, the harvested volume exceeded 5.89 million m 3 and the total annual turnover on the private forest properties was 107 million euros [9][10][11].
The forests are under the property of state and private property of individuals and legal entities (associations, schools or churches) [12].
World Wide Fund for Nature (WWF) supports the biodiversity protection of Romanian virgin forests on the UNESCO's World Heritage List [13,14]. The forests in the Carpathians represent one of the Europe's largest virgin forests, with high ecological value. Carpathian forest have been converted for, economic activities [15][16][17] and urban development [1, 11,18]. In last decades, the forests of Romania have been studied using novel methodologies, including fractal analysis [16,[19][20][21] of satellite images [11,18,22,23].
The results of these studies indicate that forests of Romanian Carpathians have suffered damage, such as diminution of forest areas and reduction of biodiversity [1, [24][25][26]. These changes resulted in a continuous decrease of forest areas and increase of urban development [27], expansion of agriculture croplands [15,28,29], climate change [15] and natural hazards [30]. Environmental sensitivity is defined in terms of various indicators that represent the physical situation of the forests (including the state property and its territorial context), in relation to the protection rules system that emphasizes its value [2].
In this study, we investigated the spatio-temporal evolution of tree cover areas, yearly deforested areas and total deforested areas, as a sum of all yearly deforested areas in the Curvature Carpathian forests in Romania by analysis of Landsat archive data from Global Forest Change (GFC) by fractal and Gray-Level Co-occurrence Matrix (GLCM) algorithms.
This report highlights the importance of forest dynamics in deforested areas and tree cover areas in Curvature Carpathians. Accordingly, our goal here was to quantify the forest diminution across the study area, from 2000 to 2014, using Landsat satellite archived data from GFC, applying image pre-processing and use of fractal and GLCM indices. Because many previous studies focused on changes in forested regions, we took a more comprehensive insight by using fractal and GLCM algorithms. The advantages of such an approach result from the compatibility of these methods to provide an added analytical value.

Study Area
The present study is focused on the Curvature Carpathians in Romania (also known as bend area of Eastern Carpathians), the last mountain division situated on the southern side of Romania. The relief spans from high mountains (Ciucas , Peak elevation: 1954 m) to depressions (Bras , ov Depression elevation: 500-600 m). The landscape changes from deep forests, in mountainous areas to agricultural land in the Bras , ov Depression. This intra-mountain depression is one of the largest in the country with Bras , ov, a city with 227,961 inhabitants in 2011 [31] (Figure 1a). Forests in Romania cover 6.86 million ha (28.8% of the national territory) and show irregular geographical distribution (57% mountain regions, 37% hills and 6% plains) [32]. Among tree species, a mosaic can be identified, with dominant broadleaved species (76%) of beech (Fagus sylvatica), fir (Abies alba), spruce (Picea abies), oak (Quercus robur), ash (Fraxinus excelsior) and pine (Pinus) [32]. In terms of deforestation, the economic need for wood, resulted in pronounced deforestation (Figure 1b). algorithms. The advantages of such an approach result from the compatibility of these methods to provide an added analytical value.

Study Area
The present study is focused on the Curvature Carpathians in Romania (also known as bend area of Eastern Carpathians), the last mountain division situated on the southern side of Romania. The relief spans from high mountains (Ciucaș Peak elevation: 1954 m) to depressions (Brașov Depression elevation: 500-600 m). The landscape changes from deep forests, in mountainous areas to agricultural land in the Brașov Depression. This intra-mountain depression is one of the largest in the country with Brașov, a city with 227,961 inhabitants in 2011 [31] (Figure 1a). Forests in Romania cover 6.86 million ha (28.8% of the national territory) and show irregular geographical distribution (57% mountain regions, 37% hills and 6% plains) [32]. Among tree species, a mosaic can be identified, with dominant broadleaved species (76%) of beech (Fagus sylvatica), fir (Abies alba), spruce (Picea abies), oak (Quercus robur), ash (Fraxinus excelsior) and pine (Pinus) [32]. In terms of deforestation, the economic need for wood, resulted in pronounced deforestation (Figure 1b).

Image Pre-Processing
GFC obtained from the Department of Geographical Sciences at University of Maryland were used for the tree cover and deforested area evaluation. This database is the result of the Landsat 7 ETM+ (30 m spatial resolution) analysis and offers the evolution of forest areas at the global level for the 2000-2014 period [33]. Images were extracted in TIFF format by use of ArcGIS from the tree cover areas and deforested areas. The initial projection of Landsat 7 ETM+ was converted from WGS 1984 World Mercator (EPSG 3395) to the Stereographic 70 coordinate system (EPSG 31700). The satellite images were transformed using ArcGIS (ESRI, California, U.S.A.) into vector (Raster to point), and the resolution of each pixel corresponds to the dimension of a point (the resolution was chosen to be 3020 × 1590 pixels-214.92 × 112.182 km). The calculation for each year from 2000 to 2014 was performed for the three types of Landsat 7 ETM+ from GFC database (forest cover loss and year of gross loss) [17]. TIFF images were binarized manually using IQM 3.5.0 (Medical University of Graz, Graz, Austria) [34] to determine the fractal indices of tree cover areas and deforested areas in ArcGIS ( Figure 2).

Image Pre-Processing
GFC obtained from the Department of Geographical Sciences at University of Maryland were used for the tree cover and deforested area evaluation. This database is the result of the Landsat 7 ETM+ (30 m spatial resolution) analysis and offers the evolution of forest areas at the global level for the 2000-2014 period [33]. Images were extracted in TIFF format by use of ArcGIS from the tree cover areas and deforested areas. The initial projection of Landsat 7 ETM+ was converted from WGS 1984 World Mercator (EPSG 3395) to the Stereographic 70 coordinate system (EPSG 31700). The satellite images were transformed using ArcGIS (ESRI, California, U.S.A.) into vector (Raster to point), and the resolution of each pixel corresponds to the dimension of a point (the resolution was chosen to be 3020 × 1590 pixels-214.92 × 112.182 km). The calculation for each year from 2000 to 2014 was performed for the three types of Landsat 7 ETM+ from GFC database (forest cover loss and year of gross loss) [17]. TIFF images were binarized manually using IQM 3.5.0 (Medical University of Graz, Graz, Austria) [34] to determine the fractal indices of tree cover areas and deforested areas in ArcGIS ( Figure 2). The changes in the forest areas, during the last decades, were analysed by fractal analysis, based on Landsat satellite images. Based on the binarized TIFF images, the evolution of fragmentation/compaction and heterogeneity/homogeneity of tree cover areas and deforested areas from Curvature Carpathians for 2000-2014 using fractal analysis for grayscale images were determined: Pyramid Dimension-PGM algorithm (PyDPGM), Fractal Fragmentation Index (FFI),
Pyramid Dimension is a fractal dimension method using image pyramids which are a sequence of identical images in different sizes. Pyramid Dimension uses pyramid images that are sequences of identical images, but at different sizes [36].
Starting from an original image I0 with the size M0 × M0 and with the gray values extended to z (x, y) ∈ [0, M0-1], an image pyramid with I0 is created as the base of the pyramid (image below). Higher The changes in the forest areas, during the last decades, were analysed by fractal analysis, based on Landsat satellite images. Based on the binarized TIFF images, the evolution of fragmentation/compaction and heterogeneity/homogeneity of tree cover areas and deforested areas from Curvature Carpathians for 2000-2014 using fractal analysis for grayscale images were determined: Pyramid Dimension-PGM algorithm (PyD PGM ), Fractal Fragmentation Index (FFI), Tug-of-War Lacunarity (Λ T−o−W ). and Cube Counting Dimension. Software IQM 3.3.0 was used for these three dimensions [34] and Gwyddion 2.4.8 (Czech Metrology Institute, Brno, Czech Republic) for Cube Counting Dimension [35]. Pyramid Dimension is a fractal dimension method using image pyramids which are a sequence of identical images in different sizes. Pyramid Dimension uses pyramid images that are sequences of identical images, but at different sizes [36]. Starting from an original image I 0 with the size M 0 × M 0 and with the gray values extended to z (x,y) ∈ [0, M 0 − 1], an image pyramid with I 0 is created as the base of the pyramid (image below). Higher images that build the I n pyramid are obtained by decreasing the size of this basic image by using different scales, ie: images In with dimensions M n × M n with M n = M 0 /S n . This downscaling can be performed with interpolation algorithms, e.g., bilinear, cubic or "neighbor nearby".

Fractal Indices Concept
In addition, downscaling can be done in two different ways. In recursive type. Each layer I n of the pyramid is chalk from the previous layer I n−1 , whereas in the base type. Each layer I n is created directly from the original image I 0 (base of the pyramid). The layers that make up the pyramid represent images investigated at multiple spatial resolutions, so they can be used to examine their fractal dimensions. One of the Pyramid Dimension approaches is Pyramid Gradient Method (PGM) and is based on area measurement approaches.
For the PGM, z is interpreted as a sampling variant of the continuous height function of the image with the shortest lateral distance h between the positions of two known values.
The area A is calculated with the equation: Different approximations can be used for the partial derivatives of Equation (1). The simplest is the approximation of the finite centered difference that can thus be calculated: where z(x, y) is the grey value of an image I with the dimension M × M at position (x, y) The derivatives must be calculated for each image point and are summed sequentially with A, to obtain the area An. In image processing, estimates for derivations are usually obtained by convolving image I with the corresponding kernels. In the analysis of PyD PGM , more advanced kernels, Sobel, are used, because it is common for the analysis of images to include the contributions of several neighboring points.
k Sobel After applying the kernels to each image of the pyramid I n the areas A n can be estimated with Equation (1). The fractal dimension PyD PGM is estimated from the slope kl of the linear distribution of log (A n ) against log (Sn) with PyD PGM = 2 − k l .
The greater the degree of space filling of an image (the image is more fragmented as the distribution of gray tones, the more uneven) the more the Pyramid Dimension fractal dimension increases and vice versa.
Cube Counting Dimension is derived directly from a definition of box-counting fractal dimension. The algorithm is based on the following steps: a cubic lattice with constant lattice l is superimposed on the z-expanded surface. Initially, l is set at X/2 (where X is the length of the edge of the surface), resulting in a lattice of 2 × 2 × 2 = 8 cubes. Then N(l) is the number of all cubes that contain at least one object pixel. The lattice constant l is then reduced stepwise by a factor of 2 and the process repeated until l equals the distance between two adjacent pixels. The slope of a plot of log(N(l)) versus log(1/l) provides the fractal dimension Df directly.
Fractal Fragmentation Index (FFI) quantifies information obtained through fractal analysis of areas and perimeters, in a single value, describing the fractal fragmentation and therefore, can be interpreted as a compaction index [20,37].
where, FFI is the fragmentation fractal index, D A is the fractal dimension of the summed areas and D P is the fractal dimension of the summed perimeters; ε represents the size of the box; logN(ε) represents the number of contiguous and non-overlapping boxes needed to cover the object area and log N (ε) represents the number of contiguous and non-overlapping boxes required for to cover only the perimeter of the object [20,37]. FFI = 0, when the analyzed fractal objects are very small, ordinarily 1-4 pixels, so that their contour cannot be extracted, D A D = D P = 0. When the areas occupied by the fractal are large and compact the value of the FFI tends to 1, and when the areas occupied by the fractal are smaller, more dispersed, more fragmented as the FFI approaches more than 0. FFI = 1 is recorded when analyzing an object Euclidean, 100% compact, without any discontinuity (D P = 1 and D A = 2).
FFI was calculated using the FFI plugin [37] according to [20]. Tug-of-War Lacunarity (Λ T−o−W ) was calculated by the frac2D plugin for IQM [38] for quantifying the degree of heterogeneity of the tree cover areas.
where N(r) = number of boxes, Z 2 = the second moment for each width as the median of s 2 values, each of which being the mean of s 1 = squares of the counter values, s 1 and s 2 are two predefined parameters that indicated the accuracy and trust.
Small values indicate the homogeneity of the spatial distribution of forests, while the large values heterogeneity.

GLCM Analysis
This is a statistical method for texture analysis that considers the spatial relationship of pixel pairs. The GLCM calculates how often a pixel with gray-level (grayscale intensity) value i occurs either horizontally, vertically, or diagonally to adjacent pixels with the value j [39]. GLCM indices were computed using ImageJ 1.5.2 software (Nationa Institute of Health, Bethesda, MD, USA), by use of the GLCM TextureToo plugin [40,41].

Calculated GLCM Features
The Angular Second Moment (Energy or Uniformity) (ASM) (Equation (6)) provides the sum of squared elements in the GLCM and measures the textural uniformity manifesting as pixel pair repetitions. It detects disorder in textures. High Energy values occur when the gray level distribution has a constant or periodic form. Energy has a normalized range. The GLCM of all homogeneous image will have many small entries.
where p(i, j) is the probability of the grey value pair i and j.
A homogeneous image will contain only a few gray levels, providing a GLCM with just a few but relatively high values of p(i, j), where i and j are coordinates of the co-occurrence matrix. Thus, the sum of squares will be high [39].
Energy values range from 0 to 1. High energy values occur when the gray level distribution has a constant or periodic shape.
Contrast (Equation (7)) measures the local variations in the GLCM (the spatial frequency of an image and is a difference moment of the GLCM). It is the difference between the highest and the lowest values of a contiguous set of pixels. It measures the number of local variations present in the image. A low contrast image presents GLCM concentration term around the principal diagonal and features low spatial frequencies.
where i and j are coordinates of the co-occurrence matrix [39]. A high contrast image indicates large differences in values between neighboring pixel sets. Correlation feature (Equation (8)) is a measure of grey tone linear dependencies in an image. It measures the joint probability occurrence of the specified pixel pair.
where i and j are coordinates of the co-occurrence matrix, whereas σ and µ represent means and standard deviations of p x and p y for a selected distance d and angle θ [39]. The correlation is higher when the image is more homogeneous. A high contrast image indicates large differences in values between neighboring pixel sets.
Inverse Difference Moment or Homogeneity (IDM) (Equation (9)) measures the closeness of the distribution of elements in the GLCM to the GLCM diagonal (measures image homogeneity as it assumes larger values for smaller gray tone differences in pair elements). It is more sensitive to the presence of near diagonal elements in the GLCM. It has a maximum value when all elements in the image have the same value. GLCM contrast and homogeneity are strongly inversely correlated in terms of equivalent distribution in the pixel pair population. It means that homogeneity decreases if contrast increases, while Energy is kept constant.
where i and j are coordinates of the co-occurrence matrix [39]. The homogeneity has a maximum value, 1, when all the elements in the image are similar (identical). Values < 1 indicate a decrease in homogeneity.
Entropy (Equation (10)) is a statistical measure of randomness and measures the disorder or complexity of an image. The entropy is large when the image is not texturally uniform and many GLCM elements have very small values. Complex textures tend to have high entropy. Entropy is strongly but inversely correlated to Energy.
where i and j are coordinates of the co-occurrence matrix [39]. The higher the entropy value, the greater the textural disorder. When the entropy is zero it means that there is no disorder.
We used in our analysis GLCM which is a second-order co-occurrence matrix because it provides information about the relative positions of the different grayscale levels in the imagery and can thus determine entropy or energy. Thus, from our point of view, it offers an advantage over the first-order occurrence matrix that provides only information about the gray-scale distribution of the image, or the Gabor Wavelet-Based texture [42] which helps to extract the characteristics of the analyzed objects (being very useful in edge, corner or blob detection) or to the Fourier band textural ordination [43,44] which allows only the extraction and characterization of the structural and textural properties such as structural heterogeneity.
Since the results of the fractal and the GLCM analyses have different units of measure and implicitly a different scale, it is very difficult to represent them on a common graph. To circumvent this issue, values were standardized according to Equation (10): where Vnom is the nominal value, Vmax is the maximum value and Vmin is the minimum value. Following the standardization of values, all data are between 0 and 1. Thus, the results are easier to plot and visualize together.

Analysis of Tree Cover Areas and Deforested Areas Using Fractal Indices
The use of fractal analysis including 2D Pyramid Dimension, 3D Cube Counting, FFI and Λ T−o−W allowed us to quantify the degree of complexity of tree cover areas and deforested areas.  Figure 3 shows that as the degree of textural fractal complexity of the tree cover areas increases, fractal fragmentation decreases. The advantage of using FFI is its information about the fragmentation degree/compaction of tree cover areas. Lower values of FFI < 0.21 indicate that Curvature Carpathians present a greater fragmentation of tree cover areas. At the level of the evolution of the tree cover areas, Λ T−o−W underlined the local effects of deforestation on the tree cover areas compaction. Due to deforestation, the lacunarity of the tree cover areas has increased, Λ T−o−W from 0.123 in 2000 to 0.106 in 2014 (Table A1) Figure 4.
The detailed analysis of Λ T−o−W of the deforested areas showed that deforestation was made in less heterogeneous areas only where the tree cover canopy is over 60%. Where the canopy is below 40%, the deforestation was more heterogeneous, where the wooded areas themselves are heterogeneous and are mainly located in depressions ( Figure 5).

GcpeGaaiiOaaaa@47EB@
</annotation> </semantics> </math> <!--MathType@End@5@5@ --> of the deforested areas showed that deforestation was made in less heterogeneous areas only where the tree cover canopy is over 60%. Where the canopy is below 40%, the deforestation was more heterogeneous, where the wooded areas themselves are heterogeneous and are mainly located in depressions ( Figure 5).  (Table A1). GGGcaaaa@46C7@ </annotation> </semantics> </math> <!--MathType@End@5@5@ --> . of regenerated areas was smaller than that of deforested areas, except where the tree cover canopy was between 60-100% (where the forest was more compact, the regeneration was lower than defrosting and made less compact) (Figure 7). Until 2004, a tendency of increasing clustering, from 0.002 to 0.01 was observed. FFI of the regenerated areas for the period 2001-2014 was very small, only 0.004 higher than the deforested areas. This fact indicates that the regeneration was accomplished more compactly than the deforestation, even if the deforestation was lower. Λ T−o−W of regenerated areas was smaller than that of deforested areas, except where the tree cover canopy was between 60-100% (where the forest was more compact, the regeneration was lower than defrosting and made less compact) ( Figure 7).
As the extent of deforestation increased during the analyzed period, Pyramid and Cube Dimensions of the Curvature Carpathian total deforested areas increased from 2001 to 2014, indicating that the deforestation was done in a chaotically manner (Figure 7). The values of Pyramid Dimension and Cube Counting of the regenerated area for 2001-2014 were lower than of the deforested area, indicating that regeneration was rather even.

GLCM Analysis
Figure 8a-c shows that the GLCM analysis provided information like the fractal texture complexity analysis obtained through Pyramid Dimension and Cube Counting. The GLCM analysis directly correlated with fractal textural complexity analysis through Entropy and Contrast and vice versa through Correlation, IDM and ASM. Thus, as the textural fractal complexity increases, the clutter increases-entropy, contrast and IDM, but uniformity and correlation decrease. GLCM confirmed that significant changes in the matrix occurred in years with more intense, homogeneous, compact deforestation (2002, 2004, 2005 and 2008) and less significant changes in years with more heterogeneous fragmented deforestation (2001, 2003 and 2007). It was also confirmed that regenerated areas in the 2001-2014 period were relatively evenly distributed when compared to deforestation. Dimension and Cube Counting of the regenerated area for 2001-2014 were lower than of the deforested area, indicating that regeneration was rather even.

GLCM Analysis
Figure 8a-c shows that the GLCM analysis provided information like the fractal texture complexity analysis obtained through Pyramid Dimension and Cube Counting. The GLCM analysis directly correlated with fractal textural complexity analysis through Entropy and Contrast and vice versa through Correlation, IDM and ASM. Thus, as the textural fractal complexity increases, the clutter increases -entropy, contrast and IDM, but uniformity and correlation decrease. GLCM confirmed that significant changes in the matrix occurred in years with more intense, homogeneous, compact deforestation (2002, 2004, 2005 and 2008) and less significant changes in years with more heterogeneous fragmented deforestation (2001, 2003 and 2007). It was also confirmed that regenerated areas in the 2001-2014 period were relatively evenly distributed when compared to deforestation.
Correlation, which represents the linear dependence of grey tones, had slightly different results compared to the other GLCM indices. The differences were more obvious for deforested areas and total deforested areas. In the total deforested areas, the linear dependence of grey tones has fallen sharply, followed by a flat region due to clustering of deforestation. Corresponding data is found in Table A2. GLCM analysis confirmed that significant changes in the matrix occurred in years with more intense, homogeneous, compact deforestation (2002,2004,2005

Discussion and Conclusions
Mapping the forest dynamics as the main objective of this study was carried out using multiple fractal and GLCM indices. Table 1 shows the main results for the deforested areas (as deforested areas for each year), total deforested areas (as cumulative deforested areas, obtained by cumulating the deforested areas for each year) and tree cover areas, from 2001 to 2014. In 2007 and 2012 the most pronounced annual average deforestation and forest area minimums were registered in 2003 and 2014, due to changes in deforestation politics which increased the permissiveness of the forest code. It is interesting to notice that the highest values of deforested areas were recorded in 2004 and 2008 when forest dynamics was influenced by legislative modifications through the law 247/2005 [45] and significant political changes. The Romanian Forestry Policy changed by introducing regulations for logging and timber harvesting, but these have not resulted in sustainable protection of forests. Table 1. Statistical summary of deforested areas/total deforested areas/tree cover areas (ha.) from 2001 to 2014 [33]. 1 the highest value of deforested areas; 2 total deforested areas; 3 the highest value of tree cover areas.
Year  Correlation, which represents the linear dependence of grey tones, had slightly different results compared to the other GLCM indices. The differences were more obvious for deforested areas and total deforested areas. In the total deforested areas, the linear dependence of grey tones has fallen sharply, followed by a flat region due to clustering of deforestation. Corresponding data is found in Table A2. GLCM analysis confirmed that significant changes in the matrix occurred in years with more intense, homogeneous, compact deforestation (2002, 2004, 2005 and 2008) and less significant changes in years with more heterogeneous fragmented deforestations (2001, 2003 and 2007).

Discussion and Conclusions
Mapping the forest dynamics as the main objective of this study was carried out using multiple fractal and GLCM indices. Table 1 shows the main results for the deforested areas (as deforested areas for each year), total deforested areas (as cumulative deforested areas, obtained by cumulating the deforested areas for each year) and tree cover areas, from 2001 to 2014. In 2007 and 2012 the most pronounced annual average deforestation and forest area minimums were registered in 2003 and 2014, due to changes in deforestation politics which increased the permissiveness of the forest code. It is interesting to notice that the highest values of deforested areas were recorded in 2004 and 2008 when forest dynamics was influenced by legislative modifications through the law 247/2005 [45] and significant political changes. The Romanian Forestry Policy changed by introducing regulations for logging and timber harvesting, but these have not resulted in sustainable protection of forests.
The evolution of the tree cover areas of the Curvature Carpathians in the period 2000-2014 revealed a general trend of decline, caused by legal and illegal logging. This situation was a consequence of the economic and legislation changes resulting in increased deforestation [16,[19][20][21]. In this study, we used fractal and GLCM analysis to study the evolution of the degree of heterogeneity and dispersion of deforested areas, total deforested areas and tree cover areas. The obtained results indicate that: (1) The tree cover areas have diminished as deforestation increases; (2) The clustering increases by deforestation (by increasing compaction and homogenization), especially after 2005; (3) Deforestation was made less compact than regeneration, although the areas were smaller; (4) The use of Pyramid Dimension allowed us to quantify the degree of textural complexity imposed by deforestation on tree cover areas; Sustainability 2019, 11, 6927 13 of 18 (5) The use of the Cube Counting Dimension allowed us to quantify the degree of textural complexity imposed by deforestation on tree cover areas through a 3D analysis; (6) GLCM analysis confirmed the fractal analysis, indicating that the link of uniformity/non-uniformity between neighboring pixels reflects the complexity of pixel relationships obtained by Pyramid Dimension and Cube Counting; (7) 3D fractal analysis reveals an added value to classical 1D and 2D fractal analyses by enabling the analysis of the complete set of data through a single fractal value; (8) Λ T−o−W revealed the local effects of deforestation on the compaction of the tree cover areas. In this regard, we show that in years with more homogeneous, compact and more intense deforestation Λ T−o−W indicated that deforestation has accentuated the process to heterogenization of the tree cover areas. In years with more heterogeneous, fragmented deforestation Λ T−o−W indicated that deforestation had limited effects on the compaction of tree cover areas; (9) In the years with minimum deforestation, FFI should have low values, close to 0, but there are also special situations, as in 2002, when the FFI, the FFI had very low values, even the deforested areas were very large (1324.5 hectares). This was due to deforestation dominated by small and numerous petitions. So, FFI captures the deforestation" behavior" better than area analysis, fractal or classical GLCM analysis. (10) Regeneration has become more heterogeneous and fragmented than deforestation. This is the first time that Pyramid Dimension and Cube Counting are analyzed for the study of forests. GLCM has also been used in the analysis of forests in the Central Group of Eastern Carpathians in Romania [46]. Compared to the previous study, our results showed that the degree of forest and deforestation disorder in the Curvature Carpathians is higher than that of the Central Group. In terms of fractal fragmentation of forests, Curvature Carpathians has proved to be more compact than the Northern Carpathian Mountains [47], Central Carpathians [46] Apuseni Mountains [32], but more fragmented than Bucegi Group [48].
Regarding Λ T−o−W , the spatial distribution of forests in the southern group proved to be more homogeneous than in the situation of the Northern Carpathian Mountains, Central Carpathians and Apuseni Mountains [32,47].
Our study has some limitations that need to be addressed. The images used, with a spatial resolution of 30 m, allowed us to capture only a coarser image of the forest models. The use of more detailed images would solve this restriction and thus improve the accuracy of fractal and GLCM analyzes. For a correct analysis all the analyzed images must have the same resolution, and the gray tones have the same characteristics for all the images.
Our findings confirm our hypothesis that fractal analysis on Landsat images from GFC, with a 30 m spatial resolution, provides valuable quantitative information on the spatiotemporal patterns of deforestation and forest areas. Thus, through fractal analysis, we have obtained complementary information on deforestation compared to classical image analysis based on image classification, because fractal analysis is able to analyze irregular spatial structures.
In conclusion, the forest change detection at a 30 m spatial resolution from GFC has helped to get a spatial and temporal dimension of forest cover change. The fractal and GLCM indices confirmed the high rate of deforestation and fragmentation of forest ecosystems from Curvature Carpathians. The trend of deforestation follows a continuous growth and causes instability of forest components and is thus unsustainable.

Acknowledgments:
The authors want to thank Mircea-Cristian Vis , an, their former colleague from Faculty of Geography, University of Bucharest for his useful advice to improve the paper.

Conflicts of Interest:
The authors declare no conflict of interest.