Processing and Assessment of Spectrometric, Stereoscopic Imagery Collected Using a Lightweight UAV Spectral Camera for Precision Agriculture

Remote Sens. 2013, 5, 5006-5039; doi:10.3390/rs5105006 Remote Sensing ISSN 2072-4292 www.mdpi.com/journal/remotesensing Article Processing and Assessment of Spectrometric, Stereoscopic Imagery Collected Using a Lightweight UAV Spectral Camera for Precision Agriculture Eija Honkavaara 1,*, Heikki Saari 2, Jere Kaivosoja 3, Ilkka Pölönen 4, Teemu Hakala 1, Paula Litkey 1, Jussi Mäkynen 2 and Liisa Pesonen 3 1 Finnish Geodetic Institute, Geodeetinrinne 2, P.O. Box 15, FI-02431 Masala, Finland; E-Mails: [email protected] (T.H.), [email protected] (P.L.) 2 VTT Photonic Devices and Measurement Solutions, P.O. Box 1000, FI-02044 VTT, Finland; E-Mails: [email protected] (H.S.), [email protected] (J.M.) 3 MTT Agrifood Research Finland, Vakolantie 55, FI-03400 Vihti, Finland; E-Mails: [email protected] (J.K.), liisa.pe[email protected] (L.P.) 4 Department of Mathematical Information Technology, University of Jyväskylä, P.O. Box 35, FI-40014 Jyväskylä, Finland; E-Mail: [email protected] * Author to whom correspondence should be addressed: E-Mail: [email protected]; Tel: +358-40-1920-835; Fax: +358-9-2955-5211. Received: 16 August 2013; in revised form: 30 September 2013 / Accepted: 6 October 2013 / Published: 14 October 2013 Abstract: Imaging using lightweight, unmanned airborne vehicles (UAVs) is one of the most rapidly developing fields in remote sensing technology. The new, tunable, Fabry-Perot interferometer-based (FPI) spectral camera, which weighs less than 700 g, makes it possible to collect spectrometric image blocks with stereoscopic overlaps using light-weight UAV platforms. This new technology is highly relevant, because it opens up new possibilities for measuring and monitoring the environment, which is becoming increasingly important for many environmental challenges. Our objectives were to investigate the processing and use of this new type of image data in precision agriculture. We developed the entire processing chain from raw images up to georeferenced reflectance images, digital surface models and biomass estimates. The processing integrates photogrammetric and quantitative remote sensing approaches. We carried out an empirical assessment using FPI spectral imagery collected at an agricultural wheat test site in the summer of 2012. Poor weather conditions during the campaign complicated the data processing, but this is one of the challenges that are faced in operational applications. The OPEN ACCESS Remote Sens. 2013, 5 5007 results indicated that the camera performed consistently and that the data processing was consistent, as well. During the agricultural experiments, promising results were obtained for biomass estimation when the spectral data was used and when an appropriate radiometric correction was applied to the data. Our results showed that the new FPI technology has a great potential in precision agriculture and indicated many possible future research topics. Keywords: photogrammetry; radiometry; spectrometry; hyperspectral; UAV; DSM; point cloud; biomass; agriculture 1. Introduction Modern airborne imaging technology based on unmanned airborne vehicles (UAVs) offers unprecedented possibilities for measuring our environment. For many applications, UAV-based airborne methods offer the possibility for cost-efficient data collection with the desired spatial and temporal resolutions. An important advantage of UAV-based technology is that the remote sensing data can be collected even under poor imaging conditions, that is, under cloud cover, which makes it truly operational in a wide range of environmental measuring applications. We focus here on lightweight systems, which is one of the most rapidly growing fields in UAV technology. The systems are quite competitive in local area applications and especially if repetitive data collection or a rapid response is needed. An appropriate sensor is a fundamental component of a UAV imaging system. The first operational, civil, lightweight UAV imaging systems typically used commercial video cameras or customer still cameras operating in selected three wide-bandwidth bands in red, green, blue and/or near-infrared spectral regions [1–3]. The recent sensor developments tailored for operation from UAVs offer enhanced possibilities for remote sensing applications in terms of better image quality, multi-spectral, hyper-spectral and thermal imaging [4–10] and laser scanning [11–13]. One interesting new sensor is a lightweight spectral camera developed by the VTT Technical Research Center of Finland (VTT). The camera is based on a piezo-actuated, Fabry-Perot interferometer (FPI) with an adjustable air gap [6,14]. This technology makes it possible to manufacture a lightweight spectral imager that can provide flexibly selectable spectral bands in a wavelength range of 400–1,000 nm. Furthermore, because the sensor produces images in a frame format, 3D information can be extracted if the images are collected with stereoscopic overlaps. In comparison to pushbroom imaging [7,9], the advantages of frame imaging include the possibility to collect image blocks with stereoscopic overlaps and the geometric and radiometric constraints provided by the rigid rectangular image geometry and multiple overlapping images. We think that this is important in particular for UAV applications, which typically utilize images collected under dynamic, vibrating and turbulent conditions. Conventional photogrammetric and remote sensing processing methods are not directly applicable for typical, small-format UAV imagery, because they have been developed for more stable data and images with a much larger spatial extent than what can be obtained with typical UAV imaging Remote Sens. 2013, 5 5008 systems. With UAV-based, small-format frame imaging, a large number—hundreds or even thousands—of overlapping images are needed to cover the desired object area. Systems are often operated under suboptimal conditions, such as below full or partial cloud cover. Despite the challenging conditions, the images must be processed accurately so that object characteristics can be interpreted on a quantitative geometric and radiometric basis using the data. Precision agriculture is one of the potential applications for hyperspectral UAV imaging [1,2,4,6,9,15–18]. In precision agriculture, the major objectives are to enable efficient use of resources, protection of the environment and documentation of applied management treatments by applying machine guidance and site-specific seeding, fertilization and plant protection. The expectation is that UAVs might provide an efficient remote sensing tool for these tasks [18]. A review by Zhang and Kovacs [16] showed that research is needed on many topics in order to develop efficient UAV-based methods for precision agriculture. In this study, we will demonstrate the use of the FPI spectral camera in a biomass estimation process for wheat crops; biomass is one of the central biophysical parameters to be estimated in precision agriculture [17]. The objectives of this investigation were to investigate a complete processing methodology for the FPI spectral imagery, as well as to demonstrate its potential in a biomass estimation process for precision agriculture. We depict a method for FPI image data processing in Section 2. We describe the test setup used for the empirical investigation in Section 3. We present the empirical results in Section 4 and discuss them in more detail in Section 5. 2. A Method for Processing FPI Spectral Data Cubes 2.1. An FPI-Based Spectral Camera The FPI-based spectral camera developed by the VTT provides a new way to collect spectrometric image blocks. The imager is based on the use of multiple orders of the Fabry-Perot interferometer together with the different spectral sensitivities of the red, green and blue pixels of the image sensor. With this arrangement, it is possible to capture three wavelength bands with a single exposure. When the FPI is placed in front of the sensor, the spectral sensitivity of each pixel is a function of the interferometer air gap. By changing the air gap, it is possible to acquire a new set of wavelengths. With smaller air gaps, it is also possible to capture only one or two wavelengths in each image. Separate short-pass and long-pass filters are needed to cut out unwanted transmissions at unused orders of the Fabry-Perot interferometer. During a flight, a predefined sequence of air gap values is applied using the FPI camera to reconstruct the spectrum for each pixel in the image (see more details in Section 2.3.3). Often, 24 different air gap values are used, and by these means, it is possible to collect 24 freely selectable spectral bands in a single flight, while the rest of the bands (0–48) are not independent. The desired spectral bands can be selected with a spectral step of 1 nm. This technology provides spectral data cubes with a rectangular image format, but each band in the data cube exposed to a different air gap value has a slightly different position and orientation. The principles of the FPI spectral camera have been described by Saari et al. [6] and Mäkynen et al. [14]. Results presented by Nackaerts et al. [19] demonstrated the feasibility of the sensor concept. Honkavaara et al. [20] developed the first Remote Sens. 2013, 5 5009 photogrammetric processing line for the FPI spectral camera, and Pölönen et al. [21] carried out the first performance assessments in precision agriculture using the 2011 prototype. The 2012 prototype was used in this investigation. It is equipped with custom optics with a focal length of 10.9 mm and an f-number of less than 3.0. The camera has a CMOSIS CMV4000 Complementary Metal Oxide Semiconductor (CMOS) red, green and blue (RGB) image sensor with an electronic shutter; the infrared cut filter has been removed from the image sensor. The sensor has a 2,048 × 2,048 pixel resolution, a pixel size of 5.5 μm and a radiometric resolution of 12 bits. In practical applications, the sensor is used in the two-times binned mode, while only part of the sensor area is used. This provides an image size of 1,024 × 648 pixels with a pixel size of 11 μm. The field of view (FOV) is ±18° in the flight direction, ±27° in the cross-flight direction and ±31° at the format corner. Application-based filters can be used, for example 500–900, 450–700, 600–1,000 or 400–500 nm filters. The spectral resolution range is 10–40 nm at the full width at half maximum (FWHM), and it is dependent on the FPI air gap value, as well as the filter selection. Table 1 shows the differences between the 2011 prototype and the 2012 prototype. Many of the parameters were improved for the 2012 prototype. The most significant improvement was the improvement of the f-number in order to improve the signal-to-noise ratio (SNR). This was achieved by improving the lens system and allowing for greater FPI ray angles (10° in comparison to 4°). The blurring of images was reduced by changing the rolling shutter to an electronic shutter. Table 1. Spectral imager specifications for 2011 and 2012 prototypes. FOV, field of view; FWHM, the full width at half maximum; FPI, Fabry-Perot interferometer. Parameter Prototype 2011 Prototype 2012 Horizontal and vertical FOV (°) >36, >26 >50, >37 Nominal focal length (mm) 9.3 10.9 Wavelength range (nm) 400–900 400–900 Spectral resolution at FWHM (nm); depending on the selection of the FPI air gap value 9–45 10–40 Spectral step (nm): adjustable by controlling the air gap of the FPI <1 <1 f-number <7 <3 Pixel size (μm); no binning/default binning 2.2/8.8 5.5/11 Maximum spectral image size (pixels) 2,592 × 1,944 2,048 × 2,048 Spectral image size with default binning (pixels) 640 × 480 1,024 × 648 Camera dimensions (mm) 65 × 65 × 130 <80 × 92 × 150 Weight (g); including battery, GPS receiver, downwelling irradiance sensors and cabling <420 <700 The sensor setup is shown in Figure 1. The entire imaging system includes the FPI spectral camera, a 32 GB compact flash memory card, irradiance sensors for measuring downwelling and upwelling irradiance (the irradiance sensor for measuring upwelling irradiance is not operational in the current setup), a GPS receiver and a lithium polymer (LiPo) battery. The system weighs less than 700 g. With this setup, more than 1,000 data cubes, each with up to 48 bands, can be collected in the two-times binned mode within a single flight. Remote Sens. 2013, 5 5010 Figure 1. Components of the FPI imaging system: FPI spectral camera, a compact flash memory card, two irradiance sensors, a GPS receiver and a LiPo battery. 2.2. FPI Spectral Camera Data Processing The FPI spectral camera data can be processed in a similar manner as small-format frame imagery, but some sensor-specific processing is also required. Figure 2 shows a general data processing chain for FPI spectral image data. The processing chain includes data collection, FPI spectral data cube generation, image orientation, digital surface model (DSM) calculations, radiometric model calculations and output product generation. When developing the data processing chain for the FPI spectral camera, our objective has been to integrate the sensor-specific processing steps into our existing processing line based on commercially available photogrammetric and remote sensing software. Figure 2. FPI spectral camera data processing chain. DSM, digital surface model. Data collection constitutes the first phase in the imaging chain. The central parameters that need to be set are the spectral sensitivities of the bands and the integration time. The spectral sensitivities are selected based on the requirements of the application. The integration time should be selected so that the images are not overexposed in relation to the bright objects and that there are good dynamics for the dark objects. Furthermore, the flight speed and flying altitude will impact the data quality (Section 2.3.3). The pre-processing phase requires that the sensor imagery be accurately corrected radiometrically based on laboratory calibrations [7,8,14]. Correcting the spectral smile and mismatch of different spectral bands requires sensor-specific processing, which are described in more detail in Section 2.3. Remote Sens. 2013, 5 5011 The objective of geometric processing is to obtain non-distorted 3D data for a desired coordinate system. The geometric processing steps involve determining the image orientations and DSM generation. An important recent advancement in photogrammetric processing is the new generation of image matching methods for DSM measurement, which is an important advantage for the image analysis process when using UAVs [20,22–27]. The objective of radiometric processing is to provide information about the object’s reflectivity [28]. With this process, all of the disturbances related to the imaging system, atmospheric influences and view/illumination-related factors need to be compensated for. We describe our approaches to radiometric processing in Section 2.4. The outputs resulting from the above-mentioned process include georeferenced reflectance image mosaics and 3D products, such as point clouds, point clouds with reflectivity information, DSMs, object models and stereomodels for visual evaluations [20]. 2.3. FPI Spectral Data Cube Generation A crucial step in the data processing chain is the construction of radiometrically and geometrically consistent spectral data cubes. The process includes three major phases: radiometric image corrections based on laboratory calibrations, spectral smile corrections and band matching. 2.3.1. Radiometric Correction Based on Laboratory Calibration The first steps in the data processing chain involve dark signal compensation and applying the sensor calibration information to the images with the usual corrections for photon response nonuniformity (PRNU) (including the CMOS array nonuniformity and lens falloff corrections). These parameters are determined in the VTT’s calibration laboratory. Dark signal compensation is carried out using a dark image collected before the image data collection. A Bayer-matrix reconstruction is carried out to obtain three-band images. The process has been described by Mäkynen et al. [14], and it will not be emphasized in this investigation. 2.3.2. Correction of the Spectral Smile Ideally, the optics of the FPI spectral camera are designed so that light rays go through the FPI as a collimated beam for a specific pixel of the image. In order to improve the f-number, this requirement was compromised, and a maximum FPI ray angle of 10° was allowed. This will cause a shift in the central wavelength of the FPI’s spectral peak (a spectral smile effect). For the most part, the peak wavelength (λ0) is linearly dependent on the cosine of the ray angle (θ) at the image plane: λ0 = 2·d·cos(θ)/m, m = 1, 2, ..., (1) where d is the air gap between the Fabry-Perot interferometer mirrors and m is the interference order at which the FPI is used. There are different approaches to handling the spectral smile: 1. Correcting images: Our assumption is that we can calculate smile-corrected images so that the corrected spectrums can be resampled from two spectrally (with a difference in Remote Sens. 2013, 5 5012 peak wavelength preferably less than 10 nm) and temporally (with a spatial displacement less than 20 pixels) adjacent image bands. 2. Using the central areas of the images: When images are collected with a minimum of 60% forward and side overlaps and when the most nadir parts of images are used, the smile effect is less than 5 nm and can be ignored in most applications (when the FWHM is 10–40 nm). 3. Resampling a “super spectrum” for each object point of the overlapping images providing variable central wavelengths. The entire spectrum can be utilized in the applications. For typical remote sensing applications and software, approaches 1 and 2, are the most functional, because they can be carried out in a separate step during the preprocessing phase. In order to utilize the entire extent of the images, approaches 1 or 3 are required. The VTT has developed a method for correcting the spectral smile by using two spatially and spectrally adjacent bands (approach 1). A correlation-based matching using shifts in the row and column directions is first carried out on the images to accurately align two of the bands. The corrected image is then formed by interpolating the desired wavelength from the adjacent bands. This simple approach is considered sufficient, because the bands that need to be matched are spectrally and spatially close to one another. 2.3.3. Band Matching The bands of a single data cube image do not overlap perfectly, due to the imaging principle of the FPI spectral camera (Section 2.1). This makes it challenging to determine the orientation of the bands, because, in principle, all of the bands and images have a different orientation. Exposing and transferring a single band to a synchronous dynamic random-access memory (SDRAM) takes approximately 75 ms, which is the time difference between two temporally adjacent bands. Recording a single cube, for example, with 24 air gaps, takes a total of 1,800 ms. The resulting spatial shift is dependent on the flight speed of the UAV, and the shift in pixel coordinates is determined by the flying height of the system. For example, with a flight speed of 5 m·s−1 and a continuous data collection mode (not stopping during image exposure), the horizontal transition of the UAV in the flight direction between temporally adjacent bands is 0.375 m, whereas it is 9 m for the entire data cube; with a flying altitude of 150 m (ground sample distance (GSD), 15 cm), the spatial displacements are thus 2.5 and 60 pixels, respectively. The UAV is also swinging and vibrating during the time of the data cube exposure, which generates nondeterministic positional and rotational mismatches for the bands. Due to these small random movements, the a priori information has to be improved by using data-based analysis. Two different approaches for processing the FPI image data are as follows: 1. Determining the orientations of and georeferencing the individual bands separately. With this approach, there are number-of-bands (typically 20–48) image blocks that need to be processed. We used this approach in our investigation with the 2011 camera prototype, where we processed five bands [20,29]. Remote Sens. 2013, 5 5013 2. Sampling the bands of individual data cubes in relation to the geometry of a reference band. Orientations are determined for the image block with reference band images, and this orientation information is then applied to all other bands. We studied this approach for this investigation. The principle aspect of approach 2 involves using one of the bands as the reference band and matching all of the other bands to it. In order to transform a band to match the geometry of the reference band, a geometric image transformation must be carried out. The parameters of the transformation are determined by utilizing tie points that are automatically measured between the bands via image matching. There are several challenges in this matching process. A scene provides different digital numbers (DN) in different spectral bands according to the spectral responses of the different objects, which complicates the image matching process. Another challenge is that, with dynamic UAV imaging, the overlaps between the different bands can be quite small. With 3D objects, a further complication is that a simple 2D modeling of the object’s geometry is not likely to provide accurate results; with homogeneous objects, the images might not contain sufficient features for matching, and repetitive features, such as those seen in crops, can cause false matches. Our current implementation combines the above two approaches. We use a few reference bands and then match several adjacent bands to these bands. By using several reference bands, we try to improve the accuracy and robustness of the image transformations. For example, we can select reference bands for each major wavelength region (blue, green, red, near-infrared (NIR)) or we could select reference bands to maximize the spatial adjacency of the bands to be matched together. We also avoid georeferencing all layers separately, which could slow down the processing in general and be too laborious in the case of challenging objects that require manual interactions. For the first implementation, we employed feature-based matching (FBM) for the band pairs using point features extracted via the Förstner operator [30]. Finally, an image transformation is carried out using an appropriate geometric model; in our system, an affine or a projective model is used. 2.4. Radiometric Correction of Frame Image Block Data Our experience is that UAV remote sensing data often need to be collected under sub-optimal illumination conditions, such as varying degrees of cloudiness. The possibility to collect data below cloud cover and under difficult conditions is also one of the major advantages of the UAV technology. We have developed two approaches for the radiometric processing of UAV image blocks in order to produce homogeneous data for non-homogeneous input data. The first approach is an image-information-based radiometric block adjustment method [20,29]. The basic principle of the approach is to use the gray values (digital numbers (DNs)) of the radiometric tie points in the overlapping images as observations and to determine the parameters of the radiometric model indirectly via the least squares principle. Currently, we use the following model for a gray value (DN): DNjk = arel_j·(aabs·Rjk(θi,θr,φ) + babs) +brel_j (2) where Rjk(θi,θr,φ) is the bi-directional reflectance factor (BRF) of the object point, k, in image j; θi and θr are the illumination and reflected light (observation) zenith angles, φi and φr are the azimuth angles, Remote Sens. 2013, 5 5014 respectively, and φ = φr − φi is the relative azimuth angle; aabs and babs are the parameters for the empirical line model for transforming the reflectance to DN and arel_j and brel_j are relative correction parameters with respect to the reference image. The impact of view/illumination geometry on the object reflectance is an important challenge in passive imaging. Our approach is to use a model-based, multiplicative, anisotropy factor. A reflectance, Rjk(θi,θr,φ), of the object point, k, can be given as follows: Rjk(θi,θr,φ) = Rk(0,0,φ)·ρmodeled(θi,θr,φ)/ρmodeled(0,0,φ) (3) where R k(0,0,φ) is the reflectance factor at the nadir viewing geometry, ρmodeled(θi,θr,φ) is the bi-directional reflectance distribution function (BRDF) for the specific view/illumination geometry and ρmodeled(0,0,φ) is the modeled reflectance in a desired nadir viewing geometry (see, also, [31]). We use the simple BRDF model developed by Walthall et al. [32]. The final BRF with coefficients a’ and b’ for the desired reflected geometry of θr = 0 is: Rjk(θi,θr,φ) = Rk(0,0,φ)(a’θr2 + b’θr cos φ + 1) (4) In the radiometric model, the absolute reflectance transformation parameters eliminate atmospheric influences. The relative additive and multiplicative parameters eliminate differences between the images that are mainly due to illumination changes and sensor instability. The BRDF model takes care of the view/illumination geometry-related issues. The numbers for the parameters are as follows: absolute calibration: 2; relative calibration: (number of images − 1) × 2; BRDF model: 2 (if the same model is used for the entire object); nadir reflectance of radiometric tie points: number of radiometric tie points. The parameters that are used depend on the conditions during imaging. This model includes many simplifications, but it can be extended using physical parameters. For this approach, a minimum of two reflectance reference targets are needed. With the current implementation, the parameters are determined separately for each band. During the image correction phase, Rk(0,0,φ) is calculated based on Equations (2) and (4). A second alternative is to utilize the irradiance measurements collected during the flight campaign, as described by Hakala et al. [33]. In this case, the relative adjustment of the images is obtained by selecting one reference image and calculating the relative multiplicative correction factors, Cj(λ), with respect to it: Cj(λ)=Eref(λ)/Ej(λ) (5) where Ej(λ) and Eref(λ) are the spectral irradiance measurements during the acquisition of image j and reference image ref. When multiplying the DNs of image j using Cj(λ), the radiometric differences caused by changes in the irradiance are eliminated. Our analysis, presented in another article [33], showed that this model is functional under certain conditions. If the assumptions are not valid, the accuracy of the correction will be reduced. 3. Empirical Investigation 3.1. Test Area An empirical campaign was carried out at the MTT Agrifood Research Finland (MTT) agricultural test site in Vihti (60°25'21''N, 24°22'28''E) (Figure 3). The test area, a 76 m by 385 m (2.9 ha) patch of Remote Sens. 2013, 5 5021 Figure 7. Examples of image quality. Bands from the left with central peak wavelengths and FWHMs in parenthesis: individual bands 7 (535.5 nm, 24.9 nm), 16 (606.2 nm, 44.0 nm) and 29 (787.5 nm, 32.1 nm) and a three-band, un-matched band composite (7, 16, 29). We estimated the signal-to-noise ratio (SNR) as a ratio of the average signal to the standard deviation in a small image window when using a tarpaulin with a nominal reflectance of 0.3 (Figure 8). This is not an accurate estimate of the SNR, because it was also influenced by the nonuniformity of the tarpaulin; still, it can be used as an indicative value, because the relationships between bands are realistic. The SNR was typically around 80, but a reduction appeared in some of the green bands (wavelength 550–600 nm), as well as in the NIR bands (wavelength >800 nm). The behavior was as expected. In the NIR band, the quantum efficiency of the CMOS image sensor decreases as the wavelengths increases. The decreasing SNR in some green bands was a result of the low transmission at those bands, due to the edge filter, which limits the spectral bandwidth and decreases the transmission. Figure 8. Signal-to-noise ratio (SNR) calculated using a tarpaulin with a nominal reflectance of 0.3 (x-axis: wavelength (nm); y-axis: SNR). 4.2. Band Matching The assessment of the band matching results indicated that the matching was successful. The numbers of tie points were 10–806, mostly >100, and the median was 276. The standard deviations of the shift parameters were 1.5, 0.9 and 0.8 pixels for the green, red and NIR reference bands, respectively. The average shifts of the bands matched to the green, red and NIR reference bands are shown in Figure 9. They were between five and −15 pixels in the flight direction and less than three pixels in the direction perpendicular to the flight direction. The maximum shifts were up to 25 pixels in the flight direction and 43 pixels in the cross-flight direction. The results presented in Figure 9 show that in the flight direction, the measured shift values were quite close to the expected values that were calculated based on the time difference between the reference band and the matched band. In the direction 0 20 40 60 80 100 120 500 550 600 650 700 750 800 850 900 SNR Wavelength (nm) Remote Sens. 2013, 5 5022 perpendicular to the flight, the shifts are shown as a function of the time difference between the bands. The average values showed minor systematic drift, which was larger as the time difference between the bands became larger; this is likely due to the possible minor drift of the camera’s x-axis with respect to the flying direction. These results show that the band matching results are consistent with our expectations. Our conclusion is that the method developed here for band matching was appropriate for the experimental testing done in this investigation. It can be improved upon in many ways if required for future applications. Efficient, automated quality control procedures can be developed by evaluating the matching and transformation statistics and the number of successful tie points and by comparing the estimated transformation parameters to the expected values. The matching methods can be further improved, as well. Figure 9. Average shift values of the affine transformation for the whole block in the flight (left) and cross-flight direction (right). Average shifts were calculated for each reference vs. band combination. 4.3. Geometric Processing 4.3.1. Orientations We determined the orientations of the reference bands in the green, red and near-infrared spectral regions (bands 7, 16 and 29). Block statistics are provided in Table 4. The standard deviations of the unit weight were 0.5–0.7 pixels, and they were best for the near-infrared band (band 29); the precision estimates for the orientation parameters also indicated better results for band 29. The better results for band 29 might be due to the higher intensity values in the NIR band (that provided a better SNR), which could provide better quality for the automatic tie point measurements (the fields are very dark on the green and red bands). On the other hand, the root-mean-square errors (RMSEs) of the GCPs indicated the poorest performance for band 29, which was probably due to the difficulties in measuring the GCPs for this band (the likely reason for this is that the properties of the paint used for the GCPs were not ideal in the NIR spectral region). The RMSEs for the targeted GCPs were on the level of one pixel, and they were expected to be a representative estimator of the georeferencing accuracy. -15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15 Expected shift (pixels) G R NIR Measured shift (pixels) -15 -10 -5 0 5 10 15 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Time difference (s) G R NIR Measured shift (pixels) Remote Sens. 2013, 5 5023 Table 4. Statistics on the block adjustments for reference bands 7, 16 and 29: standard error of unit weight (σ0), root-mean-square error (RMSE) values for the standard deviations of the perspective center positions (X0, Y0, Z0) and image rotations (ω, φ, κ) and the RMSEs for the residuals of the GCPs; the number of GCPs (N; targeted, natural). Central peak wavelengths and FWHMs of the reference bands: 7 (535.5 nm, 24.9 nm), 16 (606.2 nm, 44.0 nm) and 29 (787.5 nm, 32.1 nm). Band σ0 RMSE Positions (m) RMSE Rotations (°) RMSE GCPs (m) N (pixels) X0 Y0 Z0 ω φ κ X Y Z 7 0.68 0.28 0.27 0.06 0.105 0.107 0.019 0.08 0.10 0.05 9,10 16 0.73 0.26 0.25 0.06 0.094 0.100 0.018 0.06 0.13 0.05 11,13 29 0.49 0.24 0.24 0.06 0.091 0.093 0.016 0.14 0.13 0.03 10,13 For the self-calibrating bundle block adjustment, we estimated the principal point of autocollimation (x0, y0) and the radial distortion parameter, k1 (the radial distortion correction (dr) for radial distance r from the image center is dr = k1r3) (Table 5). The k1 parameter values were similar in the different bands, which was consistent with our expectations. This property is favorable for the band matching process. The estimated values for the principal point of autocollimation varied in the different bands; this was likely due to the non-optimality of the block for this task. The coordinates of the principal point correlated with the coordinates of the perspective center. More detailed investigations of the camera calibration should be carried out in a laboratory using suitable imaging geometry. Table 5. Estimated coordinates for the principal point of autocollimation (x0, y0) and radial distortion (k1) and their standard deviations for reference bands 7, 16 and 29. Central peak wavelengths and FWHMs of the reference bands: 7 (535.5 nm, 24.9 nm), 16 (606.2 nm, 44.0 nm), and 29 (787.5 nm, 32.1 nm). Band Parameters Standard Deviation x0 (mm) y0 (mm) k1 (mm·mm−3) x0 (mm) y0 (mm) k1 (mm·mm−3) 7 0.270 −0.023 0.00251 0.003 0.006 1.26E-05 16 0.242 −0.091 0.00252 0.002 0.005 1.07E-05 29 0.208 −0.132 0.00251 0.002 0.005 1.06E-05 The quality of the block adjustment was in accordance with our expectations. These results also indicated that the orientations of individual bands can be determined using state-of-the-art photogrammetric methods at commercial photogrammetric workstations. 4.3.2. DSM Figure 10a (left) shows a DSM obtained via automatic image matching using band 29 of the FPI image block. The matching quality was poor for tractor tracks, which appeared either as matching failures (missing points) or as outliers (high points), but otherwise, we obtained a relatively good point cloud. We also tested how to generate DSMs using the green and red reference bands, but the matching quality was poor, likely due to the lower SNR with the dark object (the fields are dark in the green and Remote Sens. 2013, 5 5024 red bands). The reference DSM extracted using the higher spatial resolution image block collected on the same day is shown in Figure 10b (left). This DSM has a higher quality, but some minor matching failures in the tractor tracks appeared with this DSM, as well. The height RMSE of the FPI spectral camera DSM was 35 cm at the GCPs; likewise, the comparison of the FPI DSM to the reference DSM at the vegetation sample locations provided an RMSE of 35 cm. Figure 10. DSM (left), vegetation heights (middle) and a scatter plot of the vegetation heights (right) as a function of the dry biomass provided using (a) FPI spectral camera images with a point interval of 20 cm and (b) high spatial resolution UAV images with a point interval of 10 cm. The dark blue areas indicate failures in matching. Airborne laser scanning data used as the bare ground information is from the National Land Survey [37]. (a) (b) We evaluated the potential for using vegetation heights taken from the DSM during the biomass estimation. We obtained the estimate for the crop height by calculating the difference between the FPI spectral camera DSM and the DSM based on airborne laser scanning (ALS DSM), which was collected in springtime on bare ground (Figure 10a (middle)). The estimated height of the vegetation was between zero and 1.4 m; the average height was 0.74 m, and the standard deviation was 0.36 m. The linear regression between the dry biomass values and the vegetation height did not show any FPI: y = 8E-05x + 0.6024 R² = 0.019 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ,0 1,000 2,000 3,000 Dry biomass (kg/ha) fpi Vegetation height (m) CIR: y = 0.0001x + 0.5575 R² = 0.3601 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ,0 1,000 2,000 3,000 Dry biomass (kg/ha) cir Vegetation height (m) Remote Sens. 2013, 5 5025 correlation (the R2 value of regression was zero) (Figure 10a (right)). In the more accurate reference DSM (Figure 10b), the corresponding values were as follows: a vegetation height of 0.44 to 0.99 m, an average height of 0.76 m and a standard deviation of 0.12 m; the R2 value for the linear regression of dry biomass and the vegetation height was 0.36. With both DSMs, the vegetation height was slightly overestimated. It was possible to visually identify the areas without vegetation using both DSMs. Using the reference DSM, the areas with 0% fertilization could also be visually detected; this was not possible with the FPI DSM. The FPI DSM was not of sufficient quality to reveal differences in the vegetation heights, which was an expected result, due to the relatively high height deviations in the DSM. With both DSMs, it is expected that the DSM is reliable only within the area surrounded by the GCPs. The reason for the worse DSM quality when using the FPI spectral camera could be due to the lower SNR and narrower bands, as well as to the fact that the matching software and the parameters might not have been ideal for this type of data; we tested several parameter options using the NGATE software, but they did not improve the results. One possible improvement could be to use more ideal band combinations during the matching phase. Furthermore, the matching algorithm could be further tuned, and different matching methods should be studied. The quality of the DSM can be considered promising and sufficient for image mosaic generation if the matching failures can be interpolated and filtered at a sufficient level of quality. 4.3.3. Image Mosaics We calculated the orthophoto mosaics using the orientation information and the DSM. The RMSE when using the GCPs as checkpoints was less than 0.2 m for the x and y coordinates. 4.4. Radiometric Processing We conducted radiometric processing for the entire block and for strip 3. In addition to the new results for strip 3, we reproduced relevant parts of the results from a previous study [33] with the full block in order to evaluate the impact of radiometric correction on the final application. Significant radiometric differences appeared within the image mosaics, due to the extremely variable illumination conditions. This was apparent especially for the full block (Figure 11b), while strip 3 was of a more uniform quality (Figure 11a). The radiometric correction greatly improved the homogeneity of the data (Figure 11c–f). The detail of a corrected mosaic (Figure 11g) shows the good quality of the data. We used the average coefficient of variation for the radiometric tie points as the indicator of the radiometric homogeneity of the block (Figure 12) (the correction methods are described in Table 3). The coefficients of variation values for the full block (Figure 12a) were 0.14–0.18 when a radiometric correction was not used and 0.05–0.12 when a radiometric correction was applied. The best results of 0.05–0.08 were obtained with the relative block adjustment (BA: relA). For strip 3 (Figure 12b), the values were 0.06–0.08 when a radiometric block adjustment was not performed. With radiometric correction, we obtained the lowest variation coefficients, approximately 0.02, for NIR bands; the variation coefficients were better than 0.03 for the green bands and better than 0.04 for the red bands (BA: relB, BRDF). The better values for the single strip are most likely due to the fact that less variability needed to be adjusted for the single strip. The correction based on the spectral irradiance Remote Sens. 2013, 5 5026 measurement on the ground (ground) appeared to provide better homogeneity than the UAV-based correction (uav). The image-based correction provided the best homogeneity, but there were some drift effects, which appeared as a slight tendency of the block to brighten in the north-east direction. The results from the previous study [33] also showed that all of the correction approaches—UAV irradiance (uav), ground irradiance (ground) and radiometric block adjustment (BA: relA)—provided relatively similar estimates of the illumination variations. Figure 11. Image mosaics with bands 29, 7 and 16. (a) Strip 3 and (b) full block without any radiometric corrections. (c) Strip 3 with radiometric block adjustment using additive and BRDF correction (BA: relB, BRDF) and mosaics (d) with corrections using the irradiance measurement in the UAV (uav), (e) with corrections using the irradiance measurement on the ground (ground) and (f) with the radiometric block adjustment with a multiplicative correction (BA: relA). (g) Detail of the corrected mosaic (colors are optimized for the image dynamic range). The mosaics of the full blocks are based on our previous study [33]. (a) (b) (c) (d) (e) (f) Remote Sens. 2013, 5 5027 Figure 11. Cont. (g) Figure 12. Average coefficients of variation data in radiometric tie points (a) for the full block [33] and (b) for strip 3 for different correction methods: no corr: no correction; uav: corrections using the irradiance measurement in the UAV; ground: corrections using the irradiance measurement on the ground; BA: relA: radiometric block adjustment with a multiplicative correction; BA: relB, BRDF: radiometric block adjustment with additive and BRDF correction. (a) (b) Six sample spectral profiles are shown in Figure 13 for different radiometric processing cases. The reflectances were taken as the median value in a 1 m by 1 m image window. The greatest differences between the samples appeared in the reflectance values in the NIR spectral region. In all cases, the NIR reflectance of samples with low biomass values (<1,000 kg·ha−1) were clearly lower than the NIR reflectance of samples with higher biomass. With the radiometric corrections, the high (>2,500 kg·ha−1) and medium (1,500–2,000 kg·ha−1) biomass samples could be separated (Figure 13b–e), which was not possible in the data without radiometric corrections (Figure 13a). There were also some concerns with the spectral profiles. The green reflectance peaks were not distinct, especially in cases without radiometric correction (Figure 13a), which indicated that the radiometric correction was not perfectly accurate in these cases. Furthermore, the reflectance values were low in 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 500 600 700 800 900 Wavelength (nm) no corr uav ground BA: relA Coefficient of variation 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 500 600 700 800 900 Wavelength (nm) strip3, no corr strip3, BA: relB, BRDF Coefficient of variation Remote Sens. 2013, 5 5028 the green and red reflectance regions, typically lower than the dark reflectance target used in radiometric calibration (Figure 13f); this means that the reflectance at green and red spectral regions was extrapolated. Radiometric calibration is very challenging at low reflectance, because small errors in calibration will cause relatively large errors in reflectance values. We did not have reference spectrums of the vegetation samples, so it was not possible to evaluate the absolute reflectance accuracy. The reflectance spectra of the three reflectance tarpaulins indicated that the spectra measured in images fitted well with the reference spectra measured in the laboratory (Figure 13f). Figure 13. Reflectance spectrums of six wheat samples with different dry biomass values (given as kg·ha−1 for each spectrum). The reflectances were taken as the median value in a 1 m by 1 m image window. Datasets: (a) strip 3 without any radiometric corrections, (b) strip 3 with radiometric block adjustment using additive and BRDF correction (BA: relB, BRDF), (c) full block with corrections using the irradiance measurement in the UAV (uav), (d) full block with corrections using the irradiance measurement on the ground (ground) and (e) full block with the radiometric block adjustment with a multiplicative correction (BA: relA); (f) reflectance spectrums of three reference tarpaulins, P05, P20 and P30, for the case BA: relA (P05-ref, P20-ref and P30-ref are reference spectrums measured in the laboratory). (a) (b) (c) (d) (e) (f) 0 0.1 0.2 0.3 0.4 0.5 500 600 700 800 900 Wavelength (nm) no corr 272 4 2512 2028 1497 99 4 632 Reflectance 0 0.1 0.2 0.3 0.4 0.5 500 600 700 800 900 Wavelength (nm) BA: relB, BRDF 272 4 2512 2028 1497 99 4 632 Reflectance 0 0.1 0.2 0.3 0.4 0.5 500 600 700 800 900 Wavelength (nm) uav 272 4 2512 2028 1497 99 4 632 Reflectance 0 0.1 0.2 0.3 0.4 0.5 500 600 700 800 900 Wavelength (nm) ground 272 4 2512 2028 1497 99 4 632 Reflectance 0 0.1 0.2 0.3 0.4 0.5 500 600 700 800 900 Wavelength (nm) BA: relA 272 4 2512 2028 1497 99 4 632 Reflectance 0 0.05 0.1 0.15 0.2 0.25 0.3 500 600 700 800 900 Wavelength (nm) P05 P20 P30 P05-ref P20-ref P30-ref Reflectance Remote Sens. 2013, 5 5029 The results of the radiometric correction were promising. In the future, it would be of interest to develop an approach that integrates the image-based information and the external irradiance measurements. Special attention must be paid to situations in which the illumination conditions change between the flight lines from sunny to diffuse (sun behind a cloud). Furthermore, a means for correcting the shadows and topographic effects should be integrated into the method. 4.5. Biomass Estimation Using Spectrometric Information from the FPI Spectral Camera We tested the performance of the reflectance output data from the FPI spectral camera in a biomass estimation process using a KNN estimator. Furthermore, we calculated Normalized Difference Vegetation Index (NDVI) (NDVI = (NIR 815.7 − R648)/(NIR 815.7 + R648)). Figure 14 shows the biomass estimation and NDVI statistics for different radiometric processing options. Figure 14. Biomass estimation statistics for different radiometric processing options: (a) strip 3 without any radiometric corrections (no corr); (b) strip 3 with radiometric block adjustment using additive and BRDF correction (BA: relB, BRDF); (c) full block with corrections using the irradiance measurement in the UAV (uav); (d) full block with corrections using the irradiance measurement on the ground (ground); and (e) full block with the radiometric block adjustment with a multiplicative correction (BA: relA). Figures from left: biomass estimate map (kg·ha−1), a scatter plot of the measured and estimated biomass values, Normalized Difference Vegetation Index (NDVI) map and a scatter plot of the NDVI with respect to measured biomass values. In the scatter plots, the measurements were carried out in areas of a size of 1 m by 1 m. (a) (b) Remote Sens. 2013, 5 5030 Figure 14. Cont. (c) (d) (e) In the biomass estimate maps, the areas that have no plants were quite visible in all cases (leftmost plots in Figure 14). The strips with 0% fertilization could be identified in most cases. The continuous lower biomass pattern in the middle of the area in the corrected dataset corresponds to a slight downhill terrain slope, which flattens towards the north and moves the valuable nutrients there. If a radiometric correction had not been carried out, the radiometric differences caused by the changes in the illumination were clearly visible in the mosaic and distorted the biomass estimates (Figure 14a). In the case of the correction based on the irradiance measurement with the UAV (Figure 14c), some strip-related artifacts (high biomass values in strips 2 and 4) appeared; we also identified these inaccuracies in the correction parameters in our recent study [33] and concluded that these inaccuracies were likely due to some shadowing effects of the irradiance sensor in image strips 1, 3 and 5 (solvable in future campaigns). The largest differences between the estimates with a correction based on the ground irradiance measurement (Figure 14d) and the image-based correction (Figure 14e) appeared in the south-west part of the mosaic, which had lower biomass estimates, and in the north-east part of the mosaic as higher biomass estimates for the image-based method. The results for the two image-based corrections were quite similar Remote Sens. 2013, 5 5037 18. Zecha, C.W.; Link, J.; Claupein, W. Mobile sensor platforms: Categorisation and research applications in precision farming. J. Sens. Sens. Syst. 2013, 2, 51–72. 19. Nackaerts, K.; Delauré, B.; Everaerts, J.; Michiels, B.; Holmlund, C.; Mäkynen, J.; Saari, H. Evaluation of a lightweigth UAS-prototype for hyperspectral imaging. Int. Arch. Photogramm. Remote Sens. Spat. Infor. Sci. 2010, 38, 478–483. 20. Honkavaara, E.; Kaivosoja, J.; Mäkynen, J.; Pellikka, I.; Pesonen, L.; Saari, H.; Salo, H.; Hakala, T.; Markelin, L.; Rosnell, T. Hyperspectral reflectance signatures and point clouds for precision agriculture by light weight UAV imaging system. ISPRS Annal. Photogramm. Remote Sens. Spat. Inf. Sci. 2012, I-7, 353–358. 21. Pölönen, I.; Salo, H.; Saari, H.; Kaivosoja, J.; Pesonen, L.; Honkavaara, E. Biomass estimator for NIR image with a few additional spectral band images taken from light UAS. Proc. SPIE 2012, 8369, doi:10.1117/12.918551. 22. Scholten, F.; Wewel, F. Digital 3D-data acquisition with the high resolution stereo camera-airborne (HRSC-A). Int. Arch. Photogramm. Remote Sens. Spat. Infor. Sci. 2000, 33, 901–908. 23. Leberl, F.; Irschara, A.; Pock, T.; Meixner, P.; Gruber, M.; Scholz, S.; Wiechert, A. Point clouds: Lidar versus 3D vision. Photogramm. Eng. Remote Sens. 2010, 76, 1123–1134. 24. Haala, N.; Hastedt, H.; Wolf, K.; Ressl, C.; Baltrusch, S. Digital photogrammetric camera evaluation—Generation of digital elevation models. Photogramm. Fernerkund. Geoinf. 2010, 2, 99–115. 25. Hirschmüller, H. Semi-Global Matching: Motivation, Development and Applications. In Photogrammetric Week 2011; Fritsch, D., Ed.; Wichmann Verlag: Heidelberg, Germany, 2011; pp. 173–184. 26. Rosnell, T.; Honkavaara, E. Point cloud generation from aerial image data acquired by a quadrocopter type micro unmanned aerial vehicle and a digital still camera. Sensors 2012, 12, 453–480. 27. Rosnell, T.; Honkavaara, E.; Nurminen, K. On geometric processing of multi-temproal image data collected by light UAV systems. Int. Arch. Photogramm. Remote Sens. Spat. Infor. Sci. 2011, 38, 63–68. 28. Schaepman-Strub, G.; Schaepman, M.E.; Painter, T.H.; Dangel, S.; Martonchik, J.V. Reflectance quantities in optical remote sensing—Definitions and case studies. Remote Sen. Environ. 2006, 103, 27–42. 29. Honkavara, E.; Hakala, T.; Saari, H.; Markelin, L.; Mäkynen, J.; Rosnell, T. A process for radiometric correction of UAV image blocks. Photogramm. Fernerkund. Geoinfor. 2012, doi: 10.1127/1432-8364/2012/0106. 30. Förstner, W.; Gülch, E. A Fast Operator for Detection and Precise Location of Distinct Points, Corners and Centers of Circular Features. In Proceedings of Intercommission Conference on Fast Processing of Photogrammetric Data, Interlaken, Swizerland, 2–4 June 1987; pp. 281–305. 31. Beisl, U. New Method for Correction of Bidirectional Effects in Hyperspectral Images. In Proceedings of Remote Sensing for Environmental Monitoring, GIS Applications, and Geology, Toulouse, France, 17 September 2001. Remote Sens. 2013, 5 5038 32. Walthall, C.L.; Norman, J.M.; Welles, J.M.; Campbell, G.; Blad, B.L. Simple equation to approximate the bidirectional reflectance from vegetative canopies and bare soil surfaces. Appl. Opt. 1985, 24, 383–387. 33. Hakala, T.; Honkavaara, E.; Saari, H.; Mäkynen, J.; Kaivosoja, J.; Pesonen, L.; Pölönen, I. Spectral imaging from UAVs under varying illumination conditions. Int. Arch. Photogramm. Remote Sens. Spat. Infor. Sci. 2013, XL-1/W2, 189–194. 34. Kuvatekniikka Oy Patrik Raski. Available online: http://kuvatekniikka.com/ (accessed on 5 September 2013). 35. Intersil ISL29004 Datasheet. Available online: http://www.intersil.com/content/dam/ Intersil/documents/fn62/fn6221.pdf (accessed on 9 October 2013). 36. Kotsiantis, S.B. Supervised machine learning: A review of classification techniques. Informatica 2007, 31, 249–268. 37. National Land Survey of Finland Open Data License. Available on line: http://www.maanmittauslaitos.fi/en/NLS_open_data_licence_version1_20120501 (accessed on 10 October 2013). 38. Chiang, K.-W.; Tsai, M.-L.; Chu, C.-H. The development of an UAV borne direct georeferenced photogrammetric platform for ground control point free applications. Sensors 2012, 12, 9161–9180. 39. Turner, D.; Lucieer, A.; Watson, C. An automated technique for generating georectified mosaics from ultra-high resolution Unmanned Aerial Vehicle (UAV) imagery, based on Structure from Motion (SfM) point clouds. Remote Sens. 2012, 4, 1392–1410. 40. Snavely, N. Bundler: Structure from Motion (SFM) for Unordered Image Collections. Available online: phototour.cs.washington.edu/bundler/ (accessed on 12 July 2013). 41. Mathews, A.J.; Jensen, J.L.R. Visualizing and quantifying vineyard canopy LAI using an Unmanned Aerial Vehicle (UAV) collected high density structure from motion point cloud. Remote Sens. 2013, 5, 2164–2183. 42. Honkavaara, E.; Arbiol, R.; Markelin, L.; Martinez, L.; Cramer, M.; Bovet, S.; Chandelier, L.; Ilves, R.; Klonus, S.; Marshal, P.; et al. Digital airborne photogrammetry—A new tool for quantitative remote sensing?—A state-of-the-Art review on radiometric aspects of digital photogrammetric images. Remote Sens. 2009, 1, 577–605. 43. Richter, R.; Schläpfer, D. Geo-atmospheric processing of airborne imaging spectrometry data. Part 2: Atmospheric/topographic correction. Int. J. Remote Sens. 2002, 23, 2631–2649. 44. Richter, R.; Kellenberger, T.; Kaufmann, H. Comparison of topographic correction methods. Remote Sens. 2009, 1, 184–196. 45. Beisl, U.; Telaar, J.; von Schönemark, M. Atmospheric Correction, Reflectance Calibration and BRDF Correction for ADS40 Image Data. In Proceedings of the XXI ISPRS Congress, Commission VII, Beijing, China, 3–11 July 2008. 46. Chandelier, L.; Martinoty, G. Radiometric aerial triangulation for the equalization of digital aerial images and orthoimages. Photogramm. Eng. Remote Sens. 2009, 75, 193–200. 47. Collings, S.; Cacetta, P.; Campbell, N.; Wu, X. Empirical models for radiometric calibration of digital aerial frame mosaics. IEEE Trans. Geosci. and Remote Sens. 2011, 49, 2573–2588. Remote Sens. 2013, 5 5039 48. López, D.H.; García, B.F.; Piqueras, J.G.; Aöcázar, G.V. An approach to the radiometric aerotriangulation of photogrammetric images. ISPRS J. Photogramm. Remote Sens. 2011, 66, 883–893. 49. Atzberger, C. Advances in remote sensing of agriculture: Context description, existing operational monitoring systems and major information needs. Remote Sens. 2013, 5, 949–981. 50. Watts, A.C.; Ambrosia, V.G.; Hinkley, E.A. Unmanned Aircraft Systems in remote sensing and scientific research: Classification and considerations of use. Remote Sens. 2012, 4, 1671–1692. © 2013 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).