Assessment of textural differentiations in forest resources in Romania using fractal analysis

: Deforestation and forest degradation have several negative effects on the environment including a loss of species habitats, disturbance of the water cycle and reduced ability to retain CO 2 , with consequences for global warming. We investigated the evolution of forest resources from development regions in Romania affected by both deforestation and reforestation using a non-Euclidean method based on fractal analysis. We calculated four fractal dimensions of forest areas: the fractal box-counting dimension of the forest areas, the fractal box-counting dimension of the dilated forest areas, the fractal dilation dimension and the box-counting dimension of the border of the dilated forest areas. Fractal analysis revealed morpho-structural and textural differentiations of forested, deforested and reforested areas in development regions with dominant mountain relief and high hills (more forested and compact organization) in comparison to the development regions dominated by plains or low hills (less forested, more fragmented with small and isolated clusters). Our analysis used the fractal analysis that has the advantage of analyzing the entire image, rather than studying local information, thereby enabling quantiﬁcation of the uniformity, fragmentation, heterogeneity and homogeneity of forests.


Introduction
Forests have been managed intensively due to increased demand for wood and cropland [1][2][3].Sustainable forest management is imperative when facing important environmental and socio-economic changes in forest resources, as pressure on the forests' resources has an adverse effect on several components of ecosystems (water pollution, reducing the biodiversity, climate change, soil erosion) [4].Specifically, long-term changes in forest environments have direct effects on climate change, carbon cycling, and carbon-nitrogen interactions, thus underlining the need for novel management strategies to maintain and preserve the global carbon storage and ecosystem roles of the forest [5,6].The European continent hosts a variety of forest biomes, from the temperate to boreal regions, generally characterized by considerable anthropogenic influence causing environmental changes [7][8][9][10].
Romania incorporates large continuous virgin forest ecosystems in the temperate zone of Europe [11,12].In the recent decades, the deforestation rate in Romania has rapidly intensified in the Carpathians and Subcarpathians through a surge in the demand from the wood industry, timber export growth and manufacturing of wood [3,13].The natural forest cover in Romania is under threat due to the high population growth rate [14].Urbanization processes (new buildings for residential, commercial, academic and business purposes) and the conversion of rural areas into urban developed areas have caused substantial changes in land use and land cover [15].
Romania is a competitive player on the European wood market.It produces a quarter of European hardwood lumber output due to it having the lowest selling prices on the continent [16].According to the National Institute of Statistics [17], the amount of roundwood produced in 2002 was 4.3 million m 3 from coniferous forests, 4 million m 3 from broadleaved forests, 1.1 million m 3 from hardwood forests and 1.1 million m 3 from softwood forests.A considerable increase in production has happened since the turn of the new millennium with 2015 values of 6 million m 3 coniferous forest, 6.6 million m 3 broadleaved forest, 5 million m 3 hardwood forest, and 1.1 million m 3 softwood forest.The Carpathian Mountains generally have a total growing stock of about 250-300 m 3 /ha for coniferous species with moderately lower the stock rates in the Western Carpathian Mountains (about 200-250 m 3 /ha) [18].The utilization rate of forest decreased to 48% in 2005 expressed as a percentage of the annual increment, a decrease from 60% in 1990 [19].
Based on the forest management literature in Romania, sustainable forestry should consider: 1.
The preventive role of forests through the stability of slopes mitigating landslides, floods and flash floods, and catastrophic inundations of river watersheds [20][21][22][23].

2.
The ecological function of forests (providing ecosystem services such as biodiversity) [24].

3.
The socio-economic trends governed by political changes after the post-socialist period [25,26] and forest restitutions [27,28] to former private owners and legal entities [29].
The assessment of forest cover change is a first necessary step to identify critical areas of land use changes [30].Methods of forest change mapping are mainly based on an assessment from remotely sensed imagery based on changes in spectral characteristics though time [30,31] or from indices designed for mapping forest changes such as the disturbance index [32][33][34].This index was applied to map forest cover changes in the protected areas of the Eastern Carpathians Mountains in Romania (Maramureş Mountains Nature Park, Rodna Mountains National Park, Călimani National Park) [28].
However, a spatially explicit identification and quantification of forest fragmentation, together with the spatial patterns and dynamics quantification during deforestation processes are needed.The powerful quantitative approach for this was presented by [35] based on the fractal analysis (fractal geometry).The fractal analysis provides a complementary approach to traditional per-pixel based assessment of deforestation and forest disturbance, by including a contextual dimension into the analysis.Fractal analysis in the form of the box-counting, pyramid dimension, Minkowski dimension, Tug-of-War dimension and Fractal Fragmentation Index (FFI) were recently used in the analysis of forests [36][37][38].In this study, we used a different approach to fractal analysis in relation to textural differentiations of forest resources by analysis of the fractal object dilation effects on fragmentation and heterogeneity of the forest spatial evolution.We applied this approach to forest resources from the eight development regions (DRs) in Romania to obtain information about the spatial effects and distribution of deforestation and reforestation.Through fractal analysis, we provide a quantification of the degree of fragmentation of forested areas.Ultimately, this research contributes by assisting institutional interventions to alleviate the imbalances that are the consequence of deforestation in Romania.

Study Area
This study includes eight administrative regions (development regions; DRs) in Romania, namely, North-East, South-East, South Muntenia, South-West Oltenia, West, North-West, Center and Bucharest-Ilfov DRs (Figure 1).In this study, we used a different approach to fractal analysis in relation to textural differentiations of forest resources by analysis of the fractal object dilation effects on fragmentation and heterogeneity of the forest spatial evolution.We applied this approach to forest resources from the eight development regions (DRs) in Romania to obtain information about the spatial effects and distribution of deforestation and reforestation.Through fractal analysis, we provide a quantification of the degree of fragmentation of forested areas.Ultimately, this research contributes by assisting institutional interventions to alleviate the imbalances that are the consequence of deforestation in Romania.

Study Area
This study includes eight administrative regions (development regions; DRs) in Romania, namely, North-East, South-East, South Muntenia, South-West Oltenia, West, North-West, Center and Bucharest-Ilfov DRs (Figure 1).The development regions of Romania were established in 1998 to facilitate the effective regional coordination in Romania following the European Union (EU) accession.According to the law 315/2004, development regions are not administrative-territorial units and are therefore not legal entities of the counties, thus posing challenges to forest management [39].Forest cover in Romania of 63,990 km 2 includes beech (Fagus sylvatica), fir (Abies alba), spruce (Picea abies), oak (Quercus robur), ash (Fraxinus excelsior) and pine (Pinus) [40].The largest forest area of 51.9% is in the Carpathian Mountains (formed by coniferous and deciduous species) followed by 37.2% in the hills (the Subcarpathians, the Getic Plateau, the Transylvanian Plateau, the Western Hills) and 10.9% in the plains (the West Plain and the Romanian Plain).Forest exploitations in Romania is in continuous growth in the Carpathian Mountains and in the hilly areas [20].The most deforested areas are in the The development regions of Romania were established in 1998 to facilitate the effective regional coordination in Romania following the European Union (EU) accession.According to the law 315/2004, development regions are not administrative-territorial units and are therefore not legal entities of the counties, thus posing challenges to forest management [39].Forest cover in Romania of 63,990 km 2 includes beech (Fagus sylvatica), fir (Abies alba), spruce (Picea abies), oak (Quercus robur), ash (Fraxinus excelsior) and pine (Pinus) [40].The largest forest area of 51.9% is in the Carpathian Mountains (formed by coniferous and deciduous species) followed by 37.2% in the hills (the Subcarpathians, the Getic Plateau, the Transylvanian Plateau, the Western Hills) and 10.9% in the plains (the West Plain and the Romanian Plain).Forest exploitations in Romania is in continuous growth in the Carpathian Mountains and in the hilly areas [20].The most deforested areas are in the Eastern Carpathians (counties Suceava (SV), Maramures , (MM), Harghita (HR), Covasna (CV), Vrancea (VN)) and in the Southern Carpathians (counties Arges , (AG), Vâlcea (VL), Hunedoara (HD)).
The analysis was conducted at the level of development regions because of the pronounced differences in the way forest resources have been managed between regions which in turn reflected the levels of deforestation and reforestation.

Image Preprocessing
Digital images obtained from the GFC (Global Forest Change) 2000-2012 database, the Department of Geographical Sciences at the Maryland University were used as input for the fractal analysis.This database is the result of analyzing Landsat 7 ETM+ (30 m spatial resolution) data of forest areas during the period 2000-2012 [41], representing global information on forest disturbance with the highest spatial resolution currently available.The Landsat spatial resolution provides sufficient spatial detail for a good fractal analysis and we used a time-series of satellite images representing the annual time steps from 2000 to 2012 to quantify and analyze areas of deforestation and forest reforestation.
The initial projection of the Landsat GFC data was converted from WGS 1984 World Mercator (EPSG 3395) to the Stereographic 70 coordinate system (EPSG 31700).Specific information about Landsat imagery used for the analysis mask is presented in Table 1.The satellite images were stored in uncompressed color TIFF format by using ArcGIS and the raster images were transformed to vector (Raster to point) with the resolution of each pixel being equivalent to the dimension of a point.To extract surfaces for each territorial-administrative unit (UAT), a spatial join was conducted, allowing for points to be grouped for each territorial-administrative unit.The calculation of forest change areas for each UAT and for each year from 2000 to 2012 was repeated for the three types of Landsat 7 ETM+ assessments available in the GFC database (forest cover loss, forest cover gain and year of gross loss).The Landsat 7 ETM+ imagery was segmented using the Color Deconvolution vector H & E DAB plugin for the ImageJ software version 1.49t [42] and subsequently converted to binary images, using a threshold of 1-255.The binarized images only contained black pixels with 0 value (background pixels, non-fractal pixels, meaning non-forest areas) and white pixels with 255 values (foreground pixels, fractal pixels, meaning forested, deforested and reforested areas).

Fractal Analysis
We conducted a fractal analysis of the digital binarized images of the forest area evolution affected by deforestation and reforestation for eight DRs.Fractal analysis was conducted because of the ability to provide additional textural information on the effects incurred by deforestation which provides complementary information on forest disturbance as compared to traditional per-pixel forest disturbance analysis.Fractalyse 2.4.software (Research Centre ThéMA, University of Franche-Comté, Besançon, France) was employed for fractal analysis [43,44].
The box-counting dimension was used for determining the fractal dimension of the forest areas (D Surf-dil0 ), the fractal dimension of the dilated forest areas (D Surf-dil3 ) and the fractal dimension of the border of the dilated forest areas (D Bord-dil3 ).This is the most used method to estimate fractal dimensions.The image is covered by a quadratic grid and the grid size ε is varied.For each value ε, the number of squares N(ε) containing any black pixels are counted.In Fractalyse 2.4, grids of different sizes are overlaid on a two-dimensional black and white image of forest areas and the number of black cells is counted for each grid size.The black cells represent the forest areas while the white ones indicate an absence of forest areas.A double logarithmic plot of N(ε), the number of occupied cells versus 1/ε is drawn and subsequently the slope of a linear regression gives an estimate of the box dimension.
The mathematical formula used for the calculation of the box counting fractal dimension was: where D B-C is the box-counting fractal dimension, ε is the side length of the box, and N(ε) is the number of non-overlapping and contiguous boxes of side ε required to cover the area of the fractal object [40,45].
As the zero limits cannot be applied to digital images, D B-C was estimated by means of the equation: where d is the slope of the graph log (N(ε)) against log (1/ε) [45].
The box-counting algorithm cannot offer a perfect coverage of a fractal pattern, so the number of occupied boxes is approximative.Hence, the estimation of the box-counting dimension requires a sufficient number of levels of analysis [43] but also the largest possible range of box sizes [46].
Four fractal dimensions were considered in this study: the fractal box-counting dimension of the forest areas (D Surf-dil0 ), the fractal box-counting dimension of the dilated forest areas (D Surf-dil3 ), the dilation fractal dimension (D dil3 ) and the fractal box-counting dimension of the border of the dilated forest areas (D Bord-dil3 ).We analyzed D Surf-dil3, D dil3 and D Bord-dil3 to determine whether large variation in the fractal dimension was a resultof a transition from the fractal organization with a lot of free spaces (D Surf-dil0 ) to another organization with less free spaces or a more compact organization (D Surf-dil3 ) (Figure 2).provides complementary information on forest disturbance as compared to traditional per-pixel forest disturbance analysis.Fractalyse 2.4.software (Research Centre ThéMA, University of Franche-Comté, Besançon, France) was employed for fractal analysis [43,44].The box-counting dimension was used for determining the fractal dimension of the forest areas (DSurf-dil0), the fractal dimension of the dilated forest areas (DSurf-dil3) and the fractal dimension of the border of the dilated forest areas (DBord-dil3).This is the most used method to estimate fractal dimensions.The image is covered by a quadratic grid and the grid size ε is varied.For each value ε, the number of squares N(ε) containing any black pixels are counted.In Fractalyse 2.4, grids of different sizes are overlaid on a two-dimensional black and white image of forest areas and the number of black cells is counted for each grid size.The black cells represent the forest areas while the white ones indicate an absence of forest areas.A double logarithmic plot of N(ε), the number of occupied cells versus 1/ε is drawn and subsequently the slope of a linear regression gives an estimate of the box dimension.
The mathematical formula used for the calculation of the box counting fractal dimension was: where DB-C is the box-counting fractal dimension, ε is the side length of the box, and N(ε) is the number of non-overlapping and contiguous boxes of side ε required to cover the area of the fractal object [40,45].As the zero limits cannot be applied to digital images, DB-C was estimated by means of the equation: ( where d is the slope of the graph log (N(ε)) against log (1/ε) [45].
The box-counting algorithm cannot offer a perfect coverage of a fractal pattern, so the number of occupied boxes is approximative.Hence, the estimation of the box-counting dimension requires a sufficient number of levels of analysis [43] but also the largest possible range of box sizes [46].
Four fractal dimensions were considered in this study: the fractal box-counting dimension of the forest areas (DSurf-dil0), the fractal box-counting dimension of the dilated forest areas (DSurf-dil3), the dilation fractal dimension (Ddil3) and the fractal box-counting dimension of the border of the dilated forest areas (DBord-dil3).We analyzed DSurf-dil3, Ddil3 and DBord-dil3 to determine whether large variation in the fractal dimension was a resultof a transition from the fractal organization with a lot of free spaces (DSurf-dil0) to another organization with less free spaces or a more compact organization (DSurf-dil3) (Figure 2).The first step of the morphological characterization of the forest patterns involves measuring the D Surf-dil0 which describes how the forest mass is concentrated across scales on a given area.D Surf-dil0 can take any value between 0 and 2. D Surf-dil0 equal to 0 corresponds to a limited case where the pattern is made up of a single point (e.g., a single patch of forest was only 0.022 km 2 ).D Surf-dil0 < 1 corresponds to a pattern made up of unconnected elements, detached clusters (the patches of forests are highly dispersed and mostly isolated from one another, similarly to a Fournier Dust).D Surf-dil0 > 1 corresponds to a mix of connected elements forming large clusters, connected elements forming small clusters, and isolated elements.When D Surf-dil0 is close to 2, the forest pattern is quite uniformly distributed following only a single scale with a high degree of homogeneity.In this case, the forest is uniformly distributed and most elements are connected to each other.The overall shape is then mainly a single large and compact cluster.An increasing D Surf-dil0 indicates a more dense and homogenous pattern of forest areas with increasing filling of the existing areas characterized by gaps.A decreasing D Surf-dil0 indicates a growing fragmentation in multi-scale clusters and a more uneven distribution of forest areas.This uneven fragmentation of deforested or reforested areas can be imposed by human activity or natural conditions by themselves.The quality of the fractal dimension estimation is controlled by using the Pearson correlation coefficient R. Estimations were accepted when R exceeded the value of 0.99 [47].If the fit between the empirical curve and the estimated curve is low, it is be concluded that either the pattern being studied is not fractal, or that it is multi-fractal [43].
The second step of the morphological characterization of the forest pattern involves a measure of the D dil3 .The principle of dilation is based on an algorithm introduced by Minkowski and Bouligand [48] to establish the dimension of an object using the measure theory and involves surroundings of each object point by additional border points with increasing size for the each iteration step.During the first dilation, each pixel is surrounded by a border of one pixel in width.Then, the reference element was a square of 3 2 pixels in size.At the second iteration step, each pixel was surrounded by a border of two pixels in width.The structuring element is then a square of 5 2 pixels in size [44].At the third iteration step, each pixel is surrounded by a border of three pixels in width and the structuring element is then a square of 7 2 pixels in size.
In the third step, we measured the box-counting dimension of forested patterns after three dilations.With three dilations, all details referring to 1.078 km 2 areas were unclear, meaning that the smallest spaces disappeared.The pattern is fully connected in a fractal manner and the fractal dimension of the dilated forested area (D Surf-dil3 ) should be higher than forest-free areas; D Surf-dil3 should also be higher than 1.If D Surf-dil3 differs substantially from D Surf-dil0 , the way of combining the smallest free spaces with the surrounding forest pixels follows a fractal organization very different from the spatial organization of the forest pattern at larger scales.On the contrary, if D Surf-dil3 is not very different from D Surf-dil0 , the smallest free spaces do not follow a multi-scale pattern; they only influence the value of D Surf-dil0 slightly by introducing noise into the measure [49].
The fourth step of the morphological characterization of the forest pattern was to measure the box-counting dimension of its dilated borders (D Bord-dil3 ).The morphology of the borders is different from that of the areas and involved different ways to describe the geometry of forest areas.We used fractal analysis of dilated borders in order to eliminate the tiniest free spaces whose spatial forest organization is not well known and tends to be rather irregular.
Because of the forests' fragmentation and heterogeneity, many borders of forest areas tend to exert irregularity with the D Bord-dil3 greater than 1.For more compact patterns, D Surf-dil3 tends to be higher than D Bord-dil3 and vice versa, because the shapes of forested patterns are more diverse in the largest counties.This may lead to slight differences in the absolute values of the fractal dimensions, without affecting the general conclusions of the present paper.

Spatial Evolution of Forest Covered Areas
Initially, the spatial evolution of forested, deforested and reforested areas was determined using ImageJ based on the Landsat GFC data (Table 2).The forest-covered areas had generally decreased during the period 2000-2012.Only for a few DRs the observed reforestation was higher than deforestation (Center DR and S-W Oltenia DR).At a national level, D Surf-dil0 of forested areas was 1.50 which indicated an increased complexity, fragmentation and irregular arrangement by predominantly connected forested clusters.Due to the accentuated non-uniformity of the spatial distribution of the reforested and deforested areas, D Surf-dil0 was as low as 0.93 for deforested and 0.88 for reforested areas respectively, indicating primarily detached forest clusters (Table 3).We observed that the West and Central development regions (DRs) had the highest values of D Surf-dil0 , due to the greater homogeneity of their forested areas, arranged in contiguous clusters, followed by North-West, West, South-West and North-East DRs.Scarcely forested DRs such as the South-East Dobrogea, South Muntenia and the Bucharest-Ilfov exerted a reduced values of D Surf-dil0 which was lower than the national average.These DRs are counties characterized by limited areas of forest with a dominant relief of plains or low hills (Figure 3a).A clear difference was observed between the values of fractal dimension for forested areas as compared to deforested/reforested (Figure 3b). between the values of fractal dimension for forested areas as compared to deforested/reforested (Figure 3b).The forests appear as small clusters, isolated patches or are distributed along the waters.Particularly, values were smaller (1.25) for the counties characterized by plains, as Brăila (BR) and Călărași (CL) and also for Bucharest (B) county (0.92) (Figure 4a).The fractal dimension of deforestation indicates the fact that the highest values of DSurf-dil0 are found in the Central, North-West and West regions which reveals a greater volume of exploited timber and that deforestation occurred in a rather unstructured manner.The lowest values of fractal dimensions from South Muntenia, South-East Dobrogea and Bucharest-Ilfov regions were caused by the limited areal extent of forests (Figure 4b).
The reforested areas display a smaller fractal dimension because these areas are smaller, but on the other hand they cover a large part of the fragmentation and irregularity induced by the deforestation.The fractal dimension DSurf-dil0 of a reforested area is higher than the deforested area in South-West Oltenia and Center DRs where important reforestation campaigns were implemented.On the contrary, differences are observed between DSurf-dil0 of deforested areas and reforested areas in the Bucharest-Ilfov region.The forests appear as small clusters, isolated patches or are distributed along the waters.Particularly, values were smaller (1.25) for the counties characterized by plains, as Brăila (BR) and Călăras , i (CL) and also for Bucharest (B) county (0.92) (Figure 4a).The reforested areas display a smaller fractal dimension because these areas are smaller, but on the other hand they cover a large part of the fragmentation and irregularity induced by the deforestation.The fractal dimension DSurf-dil0 of a reforested area is higher than the deforested area in South-West Oltenia and Center DRs where important reforestation campaigns were implemented.On the contrary, differences are observed between DSurf-dil0 of deforested areas and reforested areas in the Bucharest-Ilfov region.

Dilated Fractal Dimension (Ddil3)
For the forested areas, Ddil3 provides interesting information about the fractal dimension obtained after dilation (Figure 5 a,b).Two situations can be observed: the usual case with counties characterized by lower landscape (plains, low hills, sparsely forested counties) where the fractal dimension Ddil3 is less than DSurf-dil0 (e.g., DSurf-dil0 − Ddil3 for Galați (GL) was 0.32) and an atypical situation for counties with rugged landscapes (mountainous, hilly, densely forested) where Ddil was less than DSurf-dil0 (e.g.Ddil3 − DSurf-dil0 for Vrancea (VN) was 0.25) (Figure 6a).The main cause was that complex compact forestry areas markedly influence the Ddil computation.The fractal dimension of deforestation indicates the fact that the highest values of D Surf-dil0 are found in the Central, North-West and West regions which reveals a greater volume of exploited timber and that deforestation occurred in a rather unstructured manner.The lowest values of fractal dimensions from South Muntenia, South-East Dobrogea and Bucharest-Ilfov regions were caused by the limited areal extent of forests (Figure 4b).
The reforested areas display a smaller fractal dimension because these areas are smaller, but on the other hand they cover a large part of the fragmentation and irregularity induced by the deforestation.The fractal dimension D Surf-dil0 of a reforested area is higher than the deforested area in South-West Oltenia and Center DRs where important reforestation campaigns were implemented.On the contrary, differences are observed between D Surf-dil0 of deforested areas and reforested areas in the Bucharest-Ilfov region.

Dilated Fractal Dimension (D dil3 )
For the forested areas, D dil3 provides interesting information about the fractal dimension obtained after dilation (Figure 5a,b).Two situations can be observed: the usual case with counties characterized by lower landscape (plains, low hills, sparsely forested counties) where the fractal dimension D dil3 is less than D Surf-dil0 (e.g., D Surf-dil0 − D dil3 for Galat , i (GL) was 0.32) and an atypical situation for counties with rugged landscapes (mountainous, hilly, densely forested) where D dil was less than D Surf-dil0 (e.g., D dil3 − D Surf-dil0 for Vrancea (VN) was 0.25) (Figure 6a).The main cause was that complex compact forestry areas markedly influence the D dil computation.
obtained after dilation (Figure 5 a,b).Two situations can be observed: the usual case with counties characterized by lower landscape (plains, low hills, sparsely forested counties) where the fractal dimension Ddil3 is less than DSurf-dil0 (e.g., DSurf-dil0 − Ddil3 for Galați (GL) was 0.32) and an atypical situation for counties with rugged landscapes (mountainous, hilly, densely forested) where Ddil was less than DSurf-dil0 (e.g.Ddil3 − DSurf-dil0 for Vrancea (VN) was 0.25) (Figure 6a).The main cause was that complex compact forestry areas markedly influence the Ddil computation.The difference was larger than in the case of deforested areas due to less reforestation of forest areas (Table 4).Data source: Data were extracted from images provided by [41].
At the level of DRs, the lowest values were recorded in South Muntenia, South East Dobrogea and Bucharest-Ilfov regions, where the reforestation was less pronounced, (while the highest values were recorded in West, North-West and South-West Oltenia DRs with significant reforestation.(Figure 6c).  the relation D dil > D Surf-dil0 disappeared and was replaced by D dil < D Surf-dil0 .D Surf-dil0 − D dil3 varies between 0.2 in Bucharest (B) and 0.6 in Vâlcea (VL).At the DR level, the lowest values were found in Bucharest-Ilfov and South-Muntenia regions with the lowest extent of deforestations and in West and South-Oltenia DRs with the most deforestation.(Figure 6b).The same situation was observed in the reforested areas, where D Surf-dil0 − D dil3 varied between 0.18 in Bucharest (B) and 0.76 in Bihor (BH).The difference was larger than in the case of deforested areas due to less reforestation of forest areas (Table 4).Data source: Data were extracted from images provided by [41].
At the level of DRs, the lowest values were recorded in South Muntenia, South East Dobrogea and Bucharest-Ilfov regions, where the reforestation was less pronounced, (while the highest values were recorded in West, North-West and South-West Oltenia DRs with significant reforestation.(Figure 6c).

Fractal Dimension of the Dilated Forest Areas (D surf-dil3 )
After applying a three-fold dilation (D dil3 ) (Figure 7a), the fractal dimension D Surf-dil3 of forest areas was increased compared to D Surf-dil0 in all situations.The South-West Oltenia, North-East, North-West, Center, and West DRs with rugged landscapes, including mountains and high hills had a lower increase of D Surf-dil3 compared to D Surf-dil0 with a percentage increase with values between 9.6% and 23.9% according to the fractal analysis (Figure 7b).

Fractal Dimension of the Dilated Forest Areas (Dsurf-dil3)
After applying a three-fold dilation (Ddil3) (Figure 7a), the fractal dimension DSurf-dil3 of forest areas was increased compared to DSurf-dil0 in all situations.The South-West Oltenia, North-East, North-West, Center, and West DRs with rugged landscapes, including mountains and high hills had a lower increase of DSurf-dil3 compared to DSurf-dil0 with a percentage increase with values between 9.6% and 23.9% according to the fractal analysis (Figure 7b).For Bucharest-Ilfov, South-Muntenia, and South-East Dobrogea DRs with landscapes characterized by plains and low hills, there was a decrease of the fractal dimension by 4.5%-7.7 %.The situation was much clearer at the county level (Figure 8a).The counties with mountainous landscapes and widespread forested areas (Harghita (HR), Suceava (SV), and Caraș-Severin (CS) counties) had lower DSurf-dil3 values compared to DSurf-dil0 (between 2.8% and 4.3%).In counties dominated by landscapes of plains/low hills with scarce forested areas (Bucharest (B), Călărași (CL), For Bucharest-Ilfov, South-Muntenia, and South-East Dobrogea DRs with landscapes characterized by plains and low hills, there was a decrease of the fractal dimension by 4.5%-7.7%.The situation was much clearer at the county level (Figure 8a).The counties with mountainous landscapes and widespread forested areas (Harghita (HR), Suceava (SV), and Caras , -Severin (CS) counties) had lower D Surf-dil3 values compared to D Surf-dil0 (between 2.8% and 4.3%).In counties dominated by landscapes of plains/low hills with scarce forested areas (Bucharest (B), Călăras , i (CL), Teleorman (TR)), more significant increases of D Surf-dil3 compared to D Surf-dil0 were observed (with values between 11.3% and 36.6%).
DR (due to higher values of DSurf-dil3 in Caraș-Severin (CS) and Timiș (TM)), because the reforestation has been conducted in small and relatively uniform patches.At the county level, there was a strong correlation between DSurf-dil0 and DSurf-dil3 (R 2 = 0.83), but this correlation was weaker for deforested areas because the reforestation was distributed more chaotically.Regions with varied landscapes, including mountains, hills and highlands such as South-West Oltenia, North-West and Center showed a lower increase of DSurf-dil3 compared to DSurf-dil0 between 30.5% and 40.2% compared to DRs with landscapes of lower plains and hills such as Bucharest-Ilfov, South-East Dobrogea, and South Muntenia, where the dilation induced increase of fractal dimensions ranged between 58.7% and 211.8%.The values of DSurf-dil3 − DSurf-dil0 for reforested areas were much higher than the DSurf-dil3 − DSurf-dil0 for the deforested areas because of the high extent of fragmentation.However, the counties with landscapes dominated by mountains and dense forests (Gorj (GJ), Suceava (SV), Maramureș (MM)) showed a lower increase of DSurf-dil3 compared to DSurf-dil0 (from 35% to 18.3%).In counties with landscapes dominated by lower plains and hills with limited areas of forest (Bucharest (B), Ilfov (IF)), a more pronounced increase of DSurf-dil3 was detected compared to DSurf-dil0 with values between 61.3% and 334.1% (Figure 8c).

Fractal Dimension of the Dilated Forest Areas Borders (DBord-dil3)
Extracting the border of three-fold dilated images dil3 allowed for the calculation of the fractal dimension DBord-dil3.DBord-dil3 provided information about the fractal shape of forested area borders and it increases when was growing the irregularity of the borders (Figure 9  It follows that as the difference between D Surf-dil3 and D Surf-dil0 increased, the fractal fragmentation also increased, especially for counties characterized by plain landscapes, because the artificial "fractalization" of space as the dilation was more pronounced. At the level of DRs, D Surf-dil3 of deforested areas increased proportionally with the degree of dilation, except for a higher increase of D Surf-dil3 in South-West Oltenia towards the North-East DRs (Table 5).This may be explained by the smaller areas of forest in Vaslui (VS), Ias , i (IS) and Botos , ani (BT), mainly arranged in strips.At the county level, there was a very strong positive correlation between D Surf-dil0 and D Surf-dil3 , with R 2 = 0.99, indicating that the deforested area was relatively small, compact and unevenly distributed.Data source: Data were extracted from images provided by [41].
Center and North-West DRs with varied landscapes, including mountains, hills and highlands had a lower increase of D Surf-dil3 compared to D Surf-dil0 (an increase of 31.1% to 32.6%).Bucharest-Ilfov, South-East Dobrogea, and South-Muntenia DRs with lower plains and hills showed a dilation induced increase of the fractal dimension with values between 47.5% and 109.8%.The mountain landscape of highly deforested counties (Suceava (SV), Vrancea (VN), Harghita (HR)) showed a lower increase of D Surf-dil3 compared to D Surf-dil0 with values between 19.6% and 34% (Figure 8b).In counties dominated by the lower plains and hills landscapes, and with a low degree of forest areas such as Bucharest (B), Constant , a (CT), and Vaslui (VS) the increases of D Surf-dil3 compared to D Surf-dil0 , were larger with values between 54.3% and 161.5%.D Surf-dil3 of reforested areas increased proportionally to the dilation, with the exception of the higher increase in the West DR compared to the North-West DR (due to higher values of D Surf-dil3 in Caras , -Severin (CS) and Timis , (TM)), because the reforestation has been conducted in small and relatively uniform patches.At the county level, there was a strong correlation between D Surf-dil0 and D Surf-dil3 (R 2 = 0.83), but this correlation was weaker for deforested areas because the reforestation was distributed more chaotically.Regions with varied landscapes, including mountains, hills and highlands such as South-West Oltenia, North-West and Center showed a lower increase of D Surf-dil3 compared to D Surf-dil0 between 30.5% and 40.2% compared to DRs with landscapes of lower plains and hills such as Bucharest-Ilfov, South-East Dobrogea, and South Muntenia, where the dilation induced increase of fractal dimensions ranged between 58.7% and 211.8%.The values of D Surf-dil3 − D Surf-dil0 for reforested areas were much higher than the D Surf-dil3 − D Surf-dil0 for the deforested areas because of the high extent of fragmentation.However, the counties with landscapes dominated by mountains and dense forests (Gorj (GJ), Suceava (SV), Maramures , (MM)) showed a lower increase of D Surf-dil3 compared to D Surf-dil0 (from 35% to 18.3%).In counties with landscapes dominated by lower plains and hills with limited areas of forest (Bucharest (B), Ilfov (IF)), a more pronounced increase of D Surf-dil3 was detected compared to D Surf-dil0 with values between 61.3% and 334.1% (Figure 8c).

Fractal Dimension of the Dilated Forest Areas Borders (D Bord-dil3 )
Extracting the border of three-fold dilated images dil3 allowed for the calculation of the fractal dimension D Bord-dil3 .D Bord-dil3 provided information about the fractal shape of forested area borders and it increases when was growing the irregularity of the borders (Figure 9 a,b).
At the county level, the smallest fractal dimensions of forested areas (indicating lower fragmentation and complexity of the border) were found for forested areas of Bucharest (B), Vrancea (VN) and Caras , -Severin (CS) (very small but compact) whereas the largest fractal dimensions were found in Vaslui (VS), Galat , i (GL), Ias , i (IS) as the forested areas were distributed as elongated strips (Figure 10a).reforested areas.

Fractal Dimension of the Dilated Forest Areas Borders (DBord-dil3)
Extracting the border of three-fold dilated images dil3 allowed for the calculation of the fractal dimension DBord-dil3.DBord-dil3 provided information about the fractal shape of forested area borders and it increases when was growing the irregularity of the borders (Figure 9  At the county level, the smallest fractal dimensions of forested areas (indicating lower fragmentation and complexity of the border) were found for forested areas of Bucharest (B), Vrancea (VN) and Caraș-Severin (CS) (very small but compact) whereas the largest fractal dimensions were found in Vaslui (VS), Galați (GL), Iași (IS) as the forested areas were distributed as elongated strips (Figure 10a).
Lowest values of DBord-dil3 were recorded in South-East Dobrogea and Bucharest-Ilfov where deforestation was sporadic and relatively compact, the highest in N-E where the fragmentation of  forest areas was intense due to deforestations and South-West Oltenia.At the county level, for Bucharest, but also for Vrancea (VN), Covasna (CV), and Caraş-Severin (CS) counties for deforested areas, the lowest DBord-dil3 indicated fragmentation and less complexity due to very small but compact forested areas (where the deforestation was of a compact shape) and the highest fractal dimensions were found in Vaslui (VS), Cluj (CJ), and Iași (IS) counties (where the deforestation was shaped as elongated strips) (Figure 10b).At the level of DRs, the lowest values of DBord-dil3 were recorded in South-Muntenia, South-East Dobrogea and Bucharest-Ilfov regions, with small and relatively homogeneous deforested areas.The highest values of DBord-dil3 were registered in DRs with widespread deforestation of more irregular shape due to the configuration topography, and with forest poaching, such as the Center, West and North-West.The lowest values of DBord-dil3 of reforested areas were calculated, similarly to deforested areas, in South Muntenia, South-East Dobrogea and Bucharest-Ilfov DRs, with low and relatively homogeneous reforestation, especially by plantation.The highest values of DBord-dil3 were registered in DRs with fragmented landscapes (e.g., Center, West and North-West DRs), where self-regeneration and plantation campaigns were imposed (Table 6).
At the county level, the lowest fractal dimension values (between 0.57 and 1.14) were calculated the borders of forest areas in the regions of plains or low hills.e.g., Bucharest (B), Teleorman (TR), Botoşani (BT), and Vaslui (VS) counties, indicating their fragmentation and lower complexity.In contrast, the highest fractal dimension values of 1.36 and 1.44 were calculated for Alba (AB), Argeș (AG), Harghita (HR), and Suceava (SV) counties (Figure 9c).Lowest values of D Bord-dil3 were recorded in South-East Dobrogea and Bucharest-Ilfov where deforestation was sporadic and relatively compact, the highest in N-E where the fragmentation of forest areas was intense due to deforestations and South-West Oltenia.At the county level, for Bucharest, but also for Vrancea (VN), Covasna (CV), and Caraş-Severin (CS) counties for deforested areas, the lowest D Bord-dil3 indicated fragmentation and less complexity due to very small but compact forested areas (where the deforestation was of a compact shape) and the highest fractal dimensions were found in Vaslui (VS), Cluj (CJ), and Ias , i (IS) counties (where the deforestation was shaped as elongated strips) (Figure 10b).
At the level of DRs, the lowest values of D Bord-dil3 were recorded in South-Muntenia, South-East Dobrogea and Bucharest-Ilfov regions, with small and relatively homogeneous deforested areas.The highest values of D Bord-dil3 were registered in DRs with widespread deforestation of more irregular shape due to the configuration topography, and with forest poaching, such as the Center, West and North-West.The lowest values of D Bord-dil3 of reforested areas were calculated, similarly to deforested areas, in South Muntenia, South-East Dobrogea and Bucharest-Ilfov DRs, with low and relatively homogeneous reforestation, especially by plantation.The highest values of D Bord-dil3 were registered in DRs with fragmented landscapes (e.g., Center, West and North-West DRs), where self-regeneration and plantation campaigns were imposed (Table 6).Data source: Data were extracted from images provided by [41].
At the county level, the lowest fractal dimension values (between 0.57 and 1.14) were calculated the borders of forest areas in the regions of plains or low hills.e.g., Bucharest (B), Teleorman (TR), Botoşani (BT), and Vaslui (VS) counties, indicating their fragmentation and lower complexity.In contrast, the highest fractal dimension values of 1.36 and 1.44 were calculated for Alba (AB), Arges , (AG), Harghita (HR), and Suceava (SV) counties (Figure 9c).

Forest Changes in Romania
Our results showed important spatial forest changes in Romanian DRs during recent decades.By using a satellite based information on forest changes and fractal analysis, we provide information describing the state and changes in forest areas in all eight DRs of Romania during the period 2000-2012.The satellite based per-pixel assessment of forested, deforested and reforested areas was complemented by fractal analysis revealing morpho-structural and textural differentiations of forest disturbance in different DRs, providing a quantification of changes in the uniformity, homogeneity/heterogeneity and compaction/fragmentation of forests.
Forest areas in Romania are characterized by a general downward trend over the period covered by the study due to an increased level of legal and illegal logging [50].This development is reflected by the economic changes during this period and is spawned by changes in legislation which encouraged deforestations [11,[36][37][38].In a longer historical perspective, forest cover in Romania has been characterized by significant changes over the last 100 years, mostly due to socio-political factors that have had significant impact on the forest management [51].The 1926 to 1948 period has been characterized by a decrease in private land ownership, whereas from 1948 to 1989, in the socialist period, the forests were entirely state-owned.In the post-socialist period, an increase in private lands is a result of the three privatization laws from 1990 to 2004 [27,52].Logging was particularly intense during the Socialist-era in the 1960s and 1970s but forest harvesting declined during the 1990s as a consequence of the collapse of the timber market and lower institutional strength and weak law enforcement, triggering an increase in illegal logging at the expense of planned logging [27,51,52].Areas of primary forest in Romania covered aproximately 4000 km 2 in 1984 but this was reduced to 2185 km 2 in 2004 [28,53] and consequently the habitats for species were also diminished [24].

Usage of Fractal Analysis for Quantification of Forest Changes
We have evaluated the forest resources [41] using the box-counting method [44], determining the fractal dimension of the forest areas (D Surf-dil0 ), the fractal dimension of the dilated forest areas (D Surf-dil3 ), the fractal dimension of the border of the dilated forest areas (D Bord-dil3 ) and dilation fractal dimension (D dil3 ).Our results add new insights to forest cover change as compared to previous studies based on fractal analysis [3,11,13,[36][37][38].The findings of this study clearly confirmed the hypothesis that fractal analysis can provide an interesting new way of quantifying the degree of fragmentation and organization of forest areas characterized by deforestation and reforestation.Here, we used fractal analysis to study the characteristics of the forestry morphology for various levels of forested, deforested and reforested areas in Romania.Fractal analysis was used because of the ability of the fractal dimension to provide a quantitative description of the forest environment.By measuring the contrast of forestry patterns across scales, this method is applicable for the evaluation and characterization of forest area types.Through fractal analysis, we obtained complementary information on deforestation and reforestation/regeneration compared to classical change detection analysis based on image classification because the fractal analysis is able to analyze irregular spatial structures.
In the fractal analysis of forested areas, we used the box-counting method to describe spatial patterns [11,[36][37][38] and the dynamics of deforestation [34].In previous studies, fractal analysis of forest areas was determined by using the Fractal Fragmentation Index [13,36] or Minkowski Dimension and Pyramid Dimension [33].Until now, the dilation method was used to evaluate the quality of the commune or urban built-up areas [44] and we clarified the importance of utilizing the dilation method in understanding the textural organization of forest areas and the effects generated by deforestation and reforestation.The main findings are summarized and discussed in the following five points: 1.
We show that fractal analysis can offer valuable textural information about the spatial organizations of forest resources.Also, the fractal analysis highlights the spatial effects of deforestation on forest resources.2.
D dil3 was shown to be a useful index that indicates the degree of compactness of forested areas.D dil3 showed higher values than D surf-dil0 only in the situation of densely forested areas and with a great degree of compactness.Otherwise, for areas less forested, deforested or reforested, D dil3 was less than D surf-dil0 .

3.
The difference between D surf-dil3 and D surf-dil0 can be utilized as an index of the forest area complexity organization.This indicates that the value of D surf-dil3 − D surf-dil0 increases for more heterogeneous, fragmented and scattered forest areas.If forest areas were compact and organized in big clusters, the difference between D surf-dil3 and D surf-dil0 tended to be lower.4.
D Bord-dil3 offered important textural information about deforestation and reforestation.The more compact the deforestation or the reforestation, the lower was the value of D Bord-dil3 .

5.
The fractal analysis highlighted morpho-structural and textural differences between forested, deforested and reforested areas grouped by DRs or counties.Differences were observed between landscapes dominated by mountains and high hills (more forested and compact organization) as compared to landscapes of lower plains or hills (less forested, more fragmented, with many small and isolated clusters).Also, obvious differences were found between morpho-structural and textural patterns of forested areas (more compact) as a function of deforested and reforested areas (more isolated points, patches randomly distributed).At the level of DRs, areas of deforestation not followed by a sufficient reforestation show a change from a compact structure of forest areas to a chaotic structure with multiple free spaces that are non-forested pixels.

Perspectives for Sustainable Forest Management
The current study shows the importance of utilizing the dilation method in understanding the textural organization of forest areas and the effects generated by deforestation and reforestation.The method, based on fractal analysis, thereby shows the advantage of revealing the effects and patterns of forest disturbance altogether.Fractal analysis of both D Surf and D bord is involved in characterizing forestry morphology.
Modelling of deforestation through fractal analysis, targeting forest textural differentiation, provides interesting complementary information to traditional per-pixel deforestation mapping.By combining such methods, it is expected that decision makers can gain an improved starting point for understanding the underlying driving forces of deforestation, which is instrumental for optimal planning of sustainable forestry and future exploitation of forest resources.Specifically, the morpho-structural and textural information on forest fragmentation and heterogeneity provided by fractal analysis can be of added value in assessing the risk of landslides, floods and flash floods, which is higher in the watersheds from DRs characterized by high deforestation rates in the Carpathians and Subcarphatians [21,22].Also, forest areas have an important ecological function in relation to ecosystem services such as biodiversity.Here, increased landscape fragmentation is known to be a major determinant for a biodiversity decline [54,55] and fractal analysis offers a powerful way of quantifying landscape or forest fragmentation and assessment of temporal changes herein from remote sensing time-series.
However, the model did not sufficiently highlight the complexity of morphology of forest areas at a local scale, such as the ratio of clusters and isolated clumps of forest.This shortcoming may be resolved by a complementary model of local fractal analysis including mass-dimension (for areas) and ruler-dimension (for the border of the principal cluster).
It must be kept in mind that the Landsat based global forest cover (GFC) maps [41] used as an input for the fractal analysis are based on definitions of forest loss.Forest loss in the GFC data is defined as the disturbance or complete removal of tree cover canopy (below 25% tree canopy cover) including any conversion of natural forests by local communities, also including deforestation for establishment of plantations.Therefore, the forest mapping should be intersected with ground knowledge for the correct interpretation of patterns of forest loss.Moreover, the Landsat based forest cover maps are based on optical remote sensing that is known to suffer from shortcomings in the mapping of forest degradation, mainly due to the lack of sensitivity of optical instruments to changes within and below the forest canopy [56].
Finally, a more detailed analysis, e.g., differentiating the type of forests (deciduous or coniferous forest, private or state forests) might help in establishing of the novel forest management which would be better suited for practical implementation.

Conclusions
Fractal analysis of forest areas provides the possibility to calculate the fractal dimension both of shapes and areas.Therefore, the fractal analysis is a tool to evaluate spatially explicit geographic phenomena in more detail and can improve our knowledge of the spatial organization of disturbances in forested areas.We analyzed changes in forest resources from Romanian development regions affected by both deforestation and reforestation by using fractal analysis.Through calculation of four forest area fractal dimensions we generated the quantitative morpho-structural and textural information about the degree of uniformity, fragmentation, heterogeneity and homogeneity complementary to traditional non-textural estimates of changes in the spatial extent of deforestation and reforestation.Therefore, we recommend the application of fractal analysis to the time-series of satellite remote sensing data as a new approach to providing additional information for the improvement of forest area management.

Figure 1 .
Figure 1.Geographical study area by outlining the development regions in Romania.

Figure 1 .
Figure 1.Geographical study area by outlining the development regions in Romania.

Figure 2 .
Figure 2. Workflow of data pre-processing, GIS analysis and fractal analysis of deforestation processes for Romania.

Figure 2 .
Figure 2. Workflow of data pre-processing, GIS analysis and fractal analysis of deforestation processes for Romania.

Figure 3 .
Figure 3. (a) DSurf-dil0 of forested, deforested and reforested areas from development regions of Romania DRs of Romania.(b) Boxplots with notches represent the median, interquartile and extreme values of DSurf-dil0 of forested, deforested and reforested areas from DRs of Romania.

Figure 3 .
Figure 3. (a) D Surf-dil0 of forested, deforested and reforested areas from development regions of Romania DRs of Romania.(b) Boxplots with notches represent the median, interquartile and extreme values of D Surf-dil0 of forested, deforested and reforested areas from DRs of Romania.

Figure 5 .Figure 4 .
Figure 5. (a) Ddil3 of forested, deforested and reforested areas from development regions of Romania.(b) Boxplots with notches represent the median, interquartile range and extreme values of Ddil3 of forested, deforested and reforested areas from Romanian DRs.This fact was obvious at the level of the following DRs: Bucharest-Ilfov, South-Muntenia with

Figure 5 .Figure 5 .
Figure 5. (a) Ddil3 of forested, deforested and reforested areas from development regions of Romania.(b) Boxplots with notches represent the median, interquartile range and extreme values of Ddil3 of forested, deforested and reforested areas from Romanian DRs.This fact was obvious at the level of the following DRs: Bucharest-Ilfov, South-Muntenia with lower landscape and less forested areas revealing Ddil < DSurf-dil0.The West and Center DRs with rugged landscapes, dominated by mountains or hills revealed Ddil > DSurf-dil0.For deforested areas, the relation Ddil > DSurf-dil0 disappeared and was replaced by Ddil < DSurf-dil0.DSurf-dil0 − Ddil3 varies between 0.2 in Bucharest (B) and 0.6 in Vâlcea (VL).At the DR level, the lowest values were found in Bucharest-Ilfov and South-Muntenia regions with the lowest extent of deforestations and in West and South-Oltenia DRs with the most deforestation.(Figure 6b).The same situation was observed in the

Figure 7 .
Figure 7. (a) DSurf-dil3 of forested, deforested and reforested areas from development regions of Romania.(b) Boxplots with notches represent the median, interquartile range and extreme values of DSurf-dil3 of forested, deforested and reforested areas from DRs of Romania.

Figure 7 .
Figure 7. (a) D Surf-dil3 of forested, deforested and reforested areas from development regions of Romania; (b) Boxplots with notches represent the median, interquartile range and extreme values of D Surf-dil3 of forested, deforested and reforested areas from DRs of Romania.

Figure 9 .
Figure 9. (a) DBord-dil3 of forested, deforested and reforested areas from development regions of Romania.(b) Boxplots with notches represent the median, interquartile range and extreme values of DBord-dil3 of forested, deforested and reforested areas from DRs of Romania.

Figure 9 .
Figure 9. (a) D Bord-dil3 of forested, deforested and reforested areas from development regions of Romania; (b) Boxplots with notches represent the median, interquartile range and extreme values of D Bord-dil3 of forested, deforested and reforested areas from DRs of Romania.

Figure 10 .
Figure 10.(a) The fractal box-counting dimension of the border of the dilated forest areas (DBord-dil3) of forested areas.(b) The fractal box-counting dimension of the border of the dilated forest areas (DBord-dil3) of deforested areas.(c) The fractal box-counting dimension of the border of the dilated forest areas (DBord-dil3) of reforested areas.

Figure 10 .
Figure 10.(a) The fractal box-counting dimension of the border of the dilated forest areas (D Bord-dil3 ) of forested areas; (b) The fractal box-counting dimension of the border of the dilated forest areas (D Bord-dil3 ) of deforested areas; (c) The fractal box-counting dimension of the border of the dilated forest areas (D Bord-dil3 ) of reforested areas.

Table 2 .
Characteristics of forest resources per development region during the period 2000-2012.

Table 4 .
Fractal dimension of forested, deforested and reforested areas (Ddil3) from development regions of Romania between 2000 and 2012

Table 4 .
Fractal dimension of forested, deforested and reforested areas (D dil3 ) from development regions of Romania between 2000 and 2012.