Close-Range Photogrammetry and Infrared Imaging for Non-Invasive Honeybee Hive Population Assessment

Close-range photogrammetry and thermographic imaging techniques are used for the acquisition of all the data needed for the non-invasive assessment of a honeybee hive population. Temperature values complemented with precise 3D geometry generated using novel close-range photogrammetric and computer vision algorithms are used for the computation of the inner beehive temperature at each point of its surface. The methodology was validated through its application to three reference beehives with different population levels. The temperatures reached by the exterior surfaces of the hives showed a direct correlation with the population level. In addition, the knowledge of the 3D reality of the hives and the position of each temperature value allowed the positioning of the bee colonies without the need to open the hives. This way, the state of honeybee hives regarding the growth of population can be estimated without disturbing its natural development.


Introduction
In recent decades, the steady reduction of honeybee (Apis mellifera L.) colonies reported in Europe [1,2] has been proven as a worrying reality worldwide.In 2010, the United Nations through the UNEP (United Nations Environment Program) included the "Global Honey Bee Colony Disorders and Other Threats to Insect Pollinators" among one of the "Early Warning of Emerging Issues".According to this UNEP issue, which was based on data from the U.S. Department of Agriculture, in 2017, U.S. honey-producing colonies were reduced to less than half the number of those in 1950.During these years, the biggest drop in honeybee population coincided with the introduction of parasitic mites in 1982 [3].There are many factors affecting bee mortality: Viruses [4][5][6][7]; Nosema ceranae [8][9][10]; Varroa destructor [6,7,11]; pesticides [12,13]; effects of acaricides [14]; loss of genetic diversity [15]; loss of habitats [2]; and quick proliferation of the newly arrived invasive species, Asian hornet (Vespa velutina) [16][17][18].
Honeybee population control is essential both for environmental and socio-economic needs as well as for the answers to scientific unknowns about fauna offered by these insects [19].Regarding the first, the services provided by bees to the natural ecosystem are vital for human societies.Bees are the predominant and most economically important group of pollinators in most geographical regions.It is estimated that honeybees are responsible for 90-95% of all pollination carried out by insects [20].According to FAO (Food and Agriculture Organization of the United Nations), 71 of the 100 crop species which provide 90% of food worldwide are bee-pollinated.In Europe alone, 84% of the 264 crop species are animal-pollinated and 4000 vegetable varieties exist thanks to pollination by bees [21].The production value of one metric ton of pollinator-dependent crop is approximately five times higher than that of crop categories that do not depend on insects [22].
It is therefore crucial to make beekeeping a more attractive hobby and a less laborious profession, in order to encourage local apiculture and pollination.It is also essential to advance techniques that allow non-invasive evaluation of colonies as this will minimize the human effect in the local ecosystem and allow predictive decision making in maintenance tasks to ensure the survival of the largest number of colonies possible.
This work proposes and tests a methodology for the performance of the automatic assessment of a honeybee hive population and the estimation of its growing state and survival capacities based on the estimation and evaluation of the inner hive temperature extracted from data acquired with different imaging sensors.The methodology consists of the processing and registration of visible and thermographic images towards the generation of 4D temperature point clouds of hives, followed by the computation and evaluation of the inner hive temperature via determination of the thermal attenuation by conduction through the hive walls.
The imaging sensors used for the objective of the paper are Red Green Blue (RGB) and thermal cameras.RGB sensors are selected for the generation of accurate 3D models, while thermal cameras are applied to the measurement of surface temperature in the bee hives.Although the latter can be used for exploitation of the geometry of the object in the image [23], RGB cameras present a higher resolution that enables the achievement of higher accuracy and resolution [24].RGB and thermal data fusion have been performed with different methodologies: From the registration image-to-image [25] to registration RGB point cloud-to-image [26].In both cases, knowledge of the inner orientation parameters of the cameras is required, which is acquired through camera calibration.Camera calibration has been widely studied for RGB cameras, with a variety of calibration fields in two and three-dimensions [27,28], but also with algorithms for simultaneous calibration of the camera while performing image orientation [29].However, in the case of a thermal camera, the geometric calibration is more complex due to the difficulty associated with these cameras when attempting to distinguish between materials that are at the same temperature.For this reason, existing alternatives for thermal camera geometric calibration are either based on temperature differences [30] or emissivity differences [31].
The paper has been structured as follows: After this introduction, Section 2 includes a detailed explanation of the materials and methods used for data acquisition and processing.Section 3 deals with an explanation of the methodology presented through its application to a real case study, based on a set of honeybee hives selected as the case study; finally, Section 4 establishes the most relevant conclusions of the approach proposed.

RGB Camera
A digital single-lens reflex (DSLR) camera was used to acquire images, as shown in Figure 1, for the reconstruction of 3D point clouds, in addition to providing visual information of the state of the hive walls.The visible camera selected for this work is a Canon EOS 700D (from Canon INC., Tokyo, Japan) equipped with a 60 mm fixed focal lens.This camera has a CMOS (Complementary Metal-Oxide Semiconductor) sensor, which size is 22.3 mm × 14.9 mm.The size of the image captured with this sensor is 18.5 MP with a pixel size of 4.3 µm.The camera was calibrated prior to acquisition to allow for correction of the distortion and perspective effects from the data collected and the 3D reconstruction in the photogrammetric process.The camera calibration working in the visible band of the electromagnetic spectrum, see Table 1, is performed through the acquisition of multiple convergent images of a geometric pattern (known as a calibration grid) with different orientations of the camera.The adjustment of the perspective rays governing the position of the camera and the image in each acquisition allows for the determination of the inner orientation parameters of the camera (focal length, format size, principal point, and radial lens distortion; in this case, the decentering of the lens distortion is considered as negligible given its low value).The camera calibration was processed in the commercial photogrammetric station Photomodeler Pro5© (by Photomodeler Technologies, Vancouver, BC, Canada) which performs the automatic detection of the targets in each image, and computes and adjusts the orientation of each image.This results in the computation of the internal calibration parameters of the camera and the standard deviation of the calibration parameters as a value to test their quality [32].Thermographic data acquisition was performed with an NEC TH9260 camera, produced by Nippon Avionics CO., see Figure 1, measuring temperatures ranging from −20 °C to 60 °C, with a thermal accuracy of 0.03 °C at 30 °C (30 Hz).Its sensor is an Uncooled Focal Plane Array (UFPA), 640 × 480 (pixels) in size, capturing radiation between a wavelength of 7 μm and 14 μm.
The main advantage of this camera for its use in the presented methodology is the availability of its geometric calibration, see Table 2, so that its inner orientation parameters are known, such as the principal point, sensor size, and lens distortion coefficients.Thermographic cameras capture radiation in the infrared range of the spectrum, in contrast to RGB (Red Green Blue) cameras that work in the visible range.For this reason, the geometric calibration of the camera was performed The camera was calibrated prior to acquisition to allow for correction of the distortion and perspective effects from the data collected and the 3D reconstruction in the photogrammetric process.The camera calibration working in the visible band of the electromagnetic spectrum, see Table 1, is performed through the acquisition of multiple convergent images of a geometric pattern (known as a calibration grid) with different orientations of the camera.The adjustment of the perspective rays governing the position of the camera and the image in each acquisition allows for the determination of the inner orientation parameters of the camera (focal length, format size, principal point, and radial lens distortion; in this case, the decentering of the lens distortion is considered as negligible given its low value).The camera calibration was processed in the commercial photogrammetric station Photomodeler Pro5© (by Photomodeler Technologies, Vancouver, BC, Canada) which performs the automatic detection of the targets in each image, and computes and adjusts the orientation of each image.This results in the computation of the internal calibration parameters of the camera and the standard deviation of the calibration parameters as a value to test their quality [32].Thermographic data acquisition was performed with an NEC TH9260 camera, produced by Nippon Avionics CO., see Figure 1, measuring temperatures ranging from −20 • C to 60 • C, with a thermal accuracy of 0.03 • C at 30 • C (30 Hz).Its sensor is an Uncooled Focal Plane Array (UFPA), 640 × 480 (pixels) in size, capturing radiation between a wavelength of 7 µm and 14 µm.
The main advantage of this camera for its use in the presented methodology is the availability of its geometric calibration, see Table 2, so that its inner orientation parameters are known, such as the principal point, sensor size, and lens distortion coefficients.Thermographic cameras capture radiation in the infrared range of the spectrum, in contrast to RGB (Red Green Blue) cameras that work in the visible range.For this reason, the geometric calibration of the camera was performed using a specially designed calibration field, presented in [31], which is based on the capability of the thermographic cameras to detect objects with different apparent temperatures even if they are at the same real temperature, due to their different emissivity values.This calibration field consists of a wooden plate with a black background (high emissivity) on which foil targets are placed (low emissivity).In this case, the calibration was also processed in the commercial photogrammetric station Photomodeler Pro5©.In addition to the thermographic camera, the acquisition of thermal data was complemented with a contact thermometer, Testo 720, from Testo SE & Co. (Lenzkirch, Germany), see Table 3, used for the measurement of the temperature of the honeybee hives in certain points to perform the emissivity correction of each thermographic image.This way, all data needed for the thermographic analysis, and for the emissivity correction of the values measured in the thermographic images, were available.

Data Acquisition
Proper data acquisition planning is important to optimize the resources available, ensuring the high quality of the images and minimizing the capture time.
Regarding the RGB sensor, the planning of data acquisition was carried out based on the classical aerial photogrammetric principles [33] but adapted to the new algorithms and Structure from Motion (SfM) strategies [34].According to this, the scene must be exhaustively analyzed, since lighting conditions are important as they determine the shooting strategy and the values of exposure, aperture, and shutter speed of the camera.Thus, images should be acquired without illumination variations, avoiding overexposed areas and ensuring sharpness.In addition, an analysis of possible occlusions must be included, due to the probable presence of obstacles that will affect the protocol of multiple image acquisition and the desired overlap percentage between adjacent images.Regarding the geometric conditions of the photogrammetric data acquisition, the objective is to establish a multiple image acquisition protocol to reconstruct the object or scene of interest, guaranteeing the highest completeness and accuracy of the resulting 3D model.For that reason, the data acquisition must be performed following some geometric constraints.In this case, among all possible trajectories for the camera to acquire images for the generation of a 3D model, a ring-trajectory around the object was selected, as shown in Figure 2, maintaining a constant distance to the object under study (depth) and a specific separation (baseline) between camera positions or stations.Regarding the depth, this should be chosen according to the image scale or the resolution desired.This theoretical definition of "scale" in digital close-range photogrammetry is related to the geometric resolution of the pixel size projected over the object under study (Ground Sample Distance, GSD).This parameter can be calculated by considering the relationship between camera-object distance, the GSD, the focal length of the sensor, and the pixel size, as shown in Equation (1).Given the format difference between the thermographic and RGB sensors, each camera-object distance must be considered independently in order to ensure a proper GSD.
where f is the focal length of the sensor; D is the camera-object distance, ps is the size of the pixel and GSD is the Ground Sample Distance.
Regarding the baseline, avoiding large separations between stations leads to a high completeness of the 3D model due to the high image similarity and the good image correspondences during the image orientation phase.However, it is also necessary to consider that baselines that are too short involve the computation of nearly parallel intersections between perspective rays, which can result in a worse accuracy of the resulting 3D point coordinates.
As shown in Figure 2, at least four stations of these "ring" acquisition protocols, denominated as "master images", are complemented with four images following a "mosaic" protocol; these latest are known as "slave images".These slave images are taken by moving the camera slightly horizontally (left and right) and vertically (up and down) in relation to the master images, keeping a high overlap (80-90%).Complementary "master images", represented as black cameras in Figure 2, can be acquired with the aim of completing the geometry of the object under study, as established by its complexity.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 5 of 15 highest completeness and accuracy of the resulting 3D model.For that reason, the data acquisition must be performed following some geometric constraints.In this case, among all possible trajectories for the camera to acquire images for the generation of a 3D model, a ring-trajectory around the object was selected, as shown in Figure 2, maintaining a constant distance to the object under study (depth) and a specific separation (baseline) between camera positions or stations.Regarding the depth, this should be chosen according to the image scale or the resolution desired.This theoretical definition of "scale" in digital close-range photogrammetry is related to the geometric resolution of the pixel size projected over the object under study (Ground Sample Distance, GSD).This parameter can be calculated by considering the relationship between camera-object distance, the GSD, the focal length of the sensor, and the pixel size, as shown in Equation (1).Given the format difference between the thermographic and RGB sensors, each camera-object distance must be considered independently in order to ensure a proper GSD. ( where is the focal length of the sensor; is the camera-object distance, is the size of the pixel and is the Ground Sample Distance.Regarding the baseline, avoiding large separations between stations leads to a high completeness of the 3D model due to the high image similarity and the good image correspondences during the image orientation phase.However, it is also necessary to consider that baselines that are too short involve the computation of nearly parallel intersections between perspective rays, which can result in a worse accuracy of the resulting 3D point coordinates. As shown in Figure 2, at least four stations of these "ring" acquisition protocols, denominated as "master images", are complemented with four images following a "mosaic" protocol; these latest are known as "slave images".These slave images are taken by moving the camera slightly horizontally (left and right) and vertically (up and down) in relation to the master images, keeping a high overlap (80-90%).Complementary "master images", represented as black cameras in Figure 2, can be acquired with the aim of completing the geometry of the object under study, as established by its complexity.Regarding the thermographic sensor, there are several external parameters that influence the thermographic measurements, such as atmospheric attenuation, non-uniform heating of the object by solar incidence, and reflections produced by the radiation emitted and reflected by the surrounding elements.The existence of these influential factors requires the establishment of a Regarding the thermographic sensor, there are several external parameters that influence the thermographic measurements, such as atmospheric attenuation, non-uniform heating of the object by solar incidence, and reflections produced by the radiation emitted and reflected by the surrounding elements.The existence of these influential factors requires the establishment of a protocol for data acquisition with the aim of minimizing their effects and the acquisition of valid products for their processing and analysis.
The first rule is the adjustment of the ambient parameters in the camera for the atmospheric correction; that is, the camera-object distance, ambient temperature, and humidity, prior to data acquisition.
The second rule is the acquisition of thermographic images with an angle of incidence of 20-25 • between the normal of the object under study and the optical axis of the camera, avoiding the measurement of radiation reflected from the operator and other surfaces.
The third rule is the preparation of the data acquisition in order to achieve an approximate steady-state heat transfer through the hive envelope conditions.For this reason, data acquisition should be conducted at night after sunset or in the morning before sunrise, in order to minimize the solar gain effect over the hive walls, with no wind and no precipitation events during the previous 48 h [35,36].
Following these three rules, thermographic images were acquired in such a way that thermal information was obtained for the totality of the surface of the object under study.

3D Point Cloud Reconstruction
The image-based modeling technique consisting of the combination of photogrammetry and computer vision algorithms allows the reconstruction of dense 3D point clouds.This procedure is preferably performed on the RGB images since their higher spatial resolution in comparison to the thermographic images allows for more accurate results.It consists of two steps: Relative image orientation and dense point matching for 3D point cloud generation.The procedure was performed using GRAPHOS software, which is an integrated Photogrammetric Suite developed by the authors, which allows dense and metric 3D point clouds from terrestrial and UAV (Unmanned Aerial Vehicle) images to be obtained by enclosing robust photogrammetric and computer vision algorithms.More information about GRAPHOS can be found in [37].
The relative image orientation consists of the resolution of the spatial and angular position of the camera and each image acquired in the 3D space.This process starts with the automatic extraction and matching of image features through a variation of SIFT [38] (Scale-Invariant Feature Transform) algorithm, called ASIFT [39] (Affine Scale-Invariant Feature Transform), which provides effectiveness against other feature detection algorithms when considerable scale and rotation variations exist between images.The homologous points detected through the ASIFT algorithm are the input for the orientation and self-calibration procedure.The relative orientation is performed in two steps, first, a pairwise orientation is executed by relating the images to each other by means of the Longuet-Higgins algorithm [40] and second, this initial and relative approximation to the solution is used to perform a global orientation adjustment between all images and a camera self-calibration by means of the collinearity equations [41,42].As a result, the spatial and angular positioning of the RGB sensor was computed.
The knowledge of the position of the RGB sensor in the acquisition of each image enables the generation of a dense point cloud of the object under study (in this case, the beehive).In particular, a dense matching method through the SURE (Photogrammetric Surface Reconstruction from Imagery) implementation [43] based on the Semi-Global Matching technique (SGM) allows the generation of a dense and scaled 3D point cloud resulting from the determination of the 3D coordinates of each pixel.Regarding computational effort, the last step is the most expensive and time-consuming.
With these steps, a dense 3D point cloud with RGB information was obtained for each case study.The thermographic mapping of the 3D point cloud was solved manually through the identification of homologous entities (interest points) between each thermographic image and the 3D point cloud.This thermographic registration was obtained through the computation of the spatial resection [44] of the thermographic images.As described in [45], the spatial resection makes use of image coordinates and their homologous object coordinates to determine the position and orientation of the thermographic camera in each image acquisition.Following this basic definition, spatial resection of a single image can be extended to include the interior orientation parameters of the imaging sensor, or it can be reduced to include positional elements only (exterior orientation).
Once the exterior orientation of each thermographic image is calculated, the thermographic texture is transmitted to the 3D point cloud, see Figure 3.The result obtained is a thermographic 4D point cloud where each point (X,Y,Z) has a temperature value (T) associated and represented by the RGB conversion to the defined color map.Maximum and minimum temperatures of 18 • C and 8 • C, respectively, were selected in order to include all the temperature values presented and optimize the representation of temperature variations within the hives.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 7 of 15 a single image can be extended to include the interior orientation parameters of the imaging sensor, or it can be reduced to include positional elements only (exterior orientation).
Once the exterior orientation of each thermographic image is calculated, the thermographic texture is transmitted to the 3D point cloud, see Figure 3.The result obtained is a thermographic 4D point cloud where each point (X,Y,Z) has a temperature value (T) associated and represented by the RGB conversion to the defined color map.Maximum and minimum temperatures of 18 °C and 8 °C, respectively, were selected in order to include all the temperature values presented and optimize the representation of temperature variations within the hives.

Estimation of Interior Honeybee Hive Temperature
There are three different sources of infrared radiation captured by the thermographic camera if we assume that the object under study is opaque, implying that the transmission of IR (InfraRed) radiation can be considered as null.These are the emissions from the object under study; the radiation reflected from the surroundings (both attenuated by the atmospheric transmission) and the emission from the atmosphere between the camera and the object, as shown in Figure 4.The emissivity correction is based on Stefan-Boltzmann's law [46,47], assuming that a real body with a surface temperature measured by a thermometer can be considered as a real body, whereas it can be considered as a black body when temperature values are measured with a thermographic camera (with no emissivity correction).This way, the emissivity correction is applied as Equation ( 2): (2) where is the emissivity value for a real body, is the temperature measured with the thermographic camera, and is the temperature measured by the contact thermometer.

Estimation of Interior Honeybee Hive Temperature
There are three different sources of infrared radiation captured by the thermographic camera if we assume that the object under study is opaque, implying that the transmission of IR (InfraRed) radiation can be considered as null.These are the emissions from the object under study; the radiation reflected from the surroundings (both attenuated by the atmospheric transmission) and the emission from the atmosphere between the camera and the object, as shown in Figure 4.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 7 of 15 a single image can be extended to include the interior orientation parameters of the imaging sensor, or it can be reduced to include positional elements only (exterior orientation).
Once the exterior orientation of each thermographic image is calculated, the thermographic texture is transmitted to the 3D point cloud, see Figure 3.The result obtained is a thermographic 4D point cloud where each point (X,Y,Z) has a temperature value (T) associated and represented by the RGB conversion to the defined color map.Maximum and minimum temperatures of 18 °C and 8 °C, respectively, were selected in order to include all the temperature values presented and optimize the representation of temperature variations within the hives.

Estimation of Interior Honeybee Hive Temperature
There are three different sources of infrared radiation captured by the thermographic camera if we assume that the object under study is opaque, implying that the transmission of IR (InfraRed) radiation can be considered as null.These are the emissions from the object under study; the radiation reflected from the surroundings (both attenuated by the atmospheric transmission) and the emission from the atmosphere between the camera and the object, as shown in Figure 4.The emissivity correction is based on Stefan-Boltzmann's law [46,47], assuming that a real body with a surface temperature measured by a thermometer can be considered as a real body, whereas it can be considered as a black body when temperature values are measured with a thermographic camera (with no emissivity correction).This way, the emissivity correction is applied as Equation ( 2): (2) where is the emissivity value for a real body, is the temperature measured with the thermographic camera, and is the temperature measured by the contact thermometer.The emissivity correction is based on Stefan-Boltzmann's law [46,47], assuming that a real body with a surface temperature measured by a thermometer can be considered as a real body, whereas it can be considered as a black body when temperature values are measured with a thermographic camera (with no emissivity correction).This way, the emissivity correction is applied as Equation ( 2): where ε rb is the emissivity value for a real body, T bb is the temperature measured with the thermographic camera, and T rb is the temperature measured by the contact thermometer.
In addition, the correction of the atmospheric effect is directly performed by the thermographic camera through the configuration of the ambient parameters (temperature and relative humidity) and the camera-object distance.
Assuming that the heat transfer through the hive wall is due to the conduction stream, while the heat transfer from the outer part of the hive wall to the thermographic camera combines convective and radiant streams, see Figure 5, which can be explained with the following equations, Equations ( 3)-( 5): where Q COND , Q CONV , and Q RAD are the conductive, convective, and radiant heat transfers (in W), respectively.U is the overall heat transfer coefficient for the beehive wall (W/m 2 K), h out is the thermal convective coefficient of the outer part of the beehive wall (W/m 2 K), A is the area of the walls (discarding the area corresponding to the top cover, floor, and legs of the hive) where the heat transfer takes place, ε is the emissivity value of the outer surface of the beehive wall, and σ is the Boltzmann constant, equal to 5.67 × 10 −8 W/m 2 K4 [47].T in is the inner temperature of the hive, T wall is the temperature of the exterior surface of the beehive wall, in this case measured in the thermographies along all the beehive surface and T out is the ambient temperature.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 8 of 15 In addition, the correction of the atmospheric effect is directly performed by the thermographic camera through the configuration of the ambient parameters (temperature and relative humidity) and the camera-object distance.
Assuming that the heat transfer through the hive wall is due to the conduction stream, while the heat transfer from the outer part of the hive wall to the thermographic camera combines convective and radiant streams, see Figure 5, which can be explained with the following equations, Equations ( 3)-( 5): (3) (4)

4
(5) where , , and are the conductive, convective, and radiant heat transfers (in W), respectively.U is the overall heat transfer coefficient for the beehive wall (W/m 2 K), is the thermal convective coefficient of the outer part of the beehive wall (W/m 2 K), is the area of the walls (discarding the area corresponding to the top cover, floor, and legs of the hive) where the heat transfer takes place, is the emissivity value of the outer surface of the beehive wall, and is the Boltzmann constant, equal to 5.67 × 10 −8 W/m 2 K4 [47].
is the inner temperature of the hive, is the temperature of the exterior surface of the beehive wall, in this case measured in the thermographies along all the beehive surface and is the ambient temperature.This total heat transfer from the exterior surface of the beehive wall to the thermographic camera is assumed to be equal to the conduction heat transfer through the hive wall, making the assumption that thermal equilibrium is obtained in the hive, implying that the interior temperature is equal to the temperature of the interior surface of the wall.Consequently, making the heat flows equal [48], the inner temperature of the beehive wall can be calculated using Equation ( 6): The thermal convection coefficient, , is calculated according to [49] assuming horizontal heat flow through the beehive walls and natural convection.The emissivity value, , for the outer surface of the beehive wall is determined based on Stefan-Boltzmann's law [46,47] and the This total heat transfer from the exterior surface of the beehive wall to the thermographic camera is assumed to be equal to the conduction heat transfer through the hive wall, making the assumption that thermal equilibrium is obtained in the hive, implying that the interior temperature is equal to the temperature of the interior surface of the wall.Consequently, making the heat flows equal [48], the inner temperature of the beehive wall can be calculated using Equation ( 6): The thermal convection coefficient, h out , is calculated according to [49] assuming horizontal heat flow through the beehive walls and natural convection.The emissivity value, ε, for the outer surface of the beehive wall is determined based on Stefan-Boltzmann's law [46,47] and the temperature values T wall and T out are measured with the thermographic camera and the contact thermometer, see Equation (2).

Geometric Preprocessing
Given that the result of the 3D reconstruction is an unordered high-density point cloud with submillimetric GSD and non-uniform spatial distribution, it is necessary to apply a preprocessing filter to normalize the 3D spatial distribution of the points.A Voxel-Grid subsampling algorithm implemented by PCL (Point Cloud Library) [50] was used to down-sample the 3D point cloud allowing the accurate evaluation of the spatial distribution of the points and considerably accelerating the processing due to a significant reduction of the data volume.These subsampling and filtering processes create a "Voxel grid" (a set of tiny 3D boxes in space) over the input cloud data.Then, in each voxel (i.e., 3D box) all the points be approximated (i.e., downsampled) to their centroid.This approach is slightly slower than approximating them to the center of the voxel, but it represents the underlying surface more accurately.The result will be a homogeneous spatial distributed point cloud that allows the quantification of surfaces directly by counting points.

Honeybee Hive Population Assessment
The AFSSET (French Agency for Environmental and Occupational Health Safety) describes in [51] the mean theoretical honeybee population per hive and by season in temperate regions, as shown in Figure 6.Mortality is extremely high when activity is resumed at the end of winter and the beginning of spring coinciding with the periods of low temperatures, shortage of food reserves, and the commencement of the use of chemicals in agriculture.During winter, the bees in the colony clump together towards the middle of the hive, surrounding the queen.At this time, they allow the temperature of the hive to drop from its normal 34 • C to around 27 • C inside the cluster to conserve energy.Bees on the outer parts of the cluster, which will usually be around 9 • C, occasionally rotate positions with the bees on the more crowded inner parts, so that all the bees can keep temperatures warm enough to survive.Once the development phase starts again, the population of the colony is multiplied and the temperature of the inner part of the hive will raise back up to about 34

Geometric Preprocessing
Given that the result of the 3D reconstruction is an unordered high-density point cloud with submillimetric GSD and non-uniform spatial distribution, it is necessary to apply a preprocessing filter to normalize the 3D spatial distribution of the points.A Voxel-Grid subsampling algorithm implemented by PCL (Point Cloud Library) [50] was used to down-sample the 3D point cloud allowing the accurate evaluation of the spatial distribution of the points and considerably accelerating the processing due to a significant reduction of the data volume.These subsampling and filtering processes create a "Voxel grid" (a set of tiny 3D boxes in space) over the input cloud data.Then, in each voxel (i.e., 3D box) all the points will be approximated (i.e., downsampled) to their centroid.This approach is slightly slower than approximating them to the center of the voxel, but it represents the underlying surface more accurately.The result will be a homogeneous spatial distributed point cloud that allows the quantification of surfaces directly by counting points.

Honeybee Hive Population Assessment
The AFSSET (French Agency for Environmental and Occupational Health Safety) describes in [51] the mean theoretical honeybee population per hive and by season in temperate regions, as shown in Figure 6.Mortality is extremely high when activity is resumed at the end of winter and the beginning of spring coinciding with the periods of low temperatures, shortage of food reserves, and the commencement of the use of chemicals in agriculture.During winter, the bees in the colony clump together towards the middle of the hive, surrounding the queen.At this time, they allow the temperature of the hive to drop from its normal 34 °C to around 27 °C inside the cluster to conserve energy.Bees on the outer parts of the cluster, which will usually be around 9 °C, occasionally rotate positions with the bees on the more crowded inner parts, so that all the bees can keep temperatures warm enough to survive.Once the development phase starts again, the population of the colony is multiplied and the temperature of the inner part of the hive will raise back up to about 34 °C.Given the facts that the survival of the colony is threatened if the temperature cannot be kept above 12 °C and the impossibility of each bee's survival below 6 °C [52], the percentage of beehive wall surface below, between, and over these temperature values is used to evaluate the population state of the colony.The segmentation of surfaces corresponding to the walls of the beehive was performed through the implementation of a pass-through filter.A pass-through filter is a simple filtering algorithm along a specified dimension cutting off values that are either inside or outside a given user range.Given the facts that the survival of the colony is threatened if the temperature cannot be kept above 12 • C and the impossibility of each bee's survival below 6 • C [52], the percentage of beehive wall surface below, between, and over these temperature values is used to evaluate the population state of the colony.The segmentation of surfaces corresponding to the walls of the beehive was performed through the implementation of a pass-through filter.A pass-through filter is a simple filtering algorithm along a specified dimension cutting off values that are either inside or outside a given user range.
In order to compare results and validate the methodology, each representative hive was opened to check their state at the moment of data acquisition.The process consists of a manual inspection, conducted by opening the top of the hive and looking between frames.Some of them will be found to be populated by bees, while others are largely empty.Due to the impossibility of accurately quantifying the number of bees in a colony, the methodology used by some authors such as [53] was applied, based on the correlation of the results with the frame count; referring to "frame count" as the number of frames full of bees, rounded to the nearest integer for each hive body.

Study Case
The proposed methodology was validated selecting three representative hives regarding their population (i.e., very low, normal, and high population with respect to the population values for each phase, shown in Figures 6-8); these hives were utilized for the production of honey for self-consumption and located in the town of Castropol (Northern Spain) (Latitude 43 • 31 N, Longitude 6 • 59 W).The first hive, shown on the left in both Figures 7 and 8, was chosen due to its low population level.The visual inspection carried out a few hours after the imagery acquisition confirmed the existence of a very small number of bees partially covering only one frame.The survival chances of colonies like this are very limited.In this case, all bees had disappeared a few days after the data acquisition.The second hive, shown in the middle image featured in Figures 7 and 8, was selected as a representative for those hives where routine monitoring is necessary.Their population level is low, partially covering less than a half of the brood chamber stacks.If their population level continues decreasing, maintenance tasks would be required to try to prevent the extinction of the colony.The third hive, shown on the right in Figures 7 and 8, was selected as a representation of strong colonies with very high populations.In this case, almost all frames are fully covered by bees producing large crowds at the top of the stacks due to the high density of bees inside the hive.In this case, the hive does not need immediate attention, being able to face the wintering phase successfully.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 10 of 15 conducted by opening the top of the hive and looking between frames.Some of them will be found to be populated by bees, while others are largely empty.Due to the impossibility of accurately quantifying the number of bees in a colony, the methodology used by some authors such as [53] was applied, based on the correlation of the results with the frame count; referring to "frame count" as the number of frames full of bees, rounded to the nearest integer for each hive body.

Study Case
The proposed methodology was validated selecting three representative hives regarding their population (i.e., very low, normal, and high population with respect to the population values for each phase, shown in Figures 6-8); these hives were utilized for the production of honey for selfconsumption and located in the town of Castropol (Northern Spain) (Latitude 43°31′ N, Longitude 6°59′ W).The first hive, shown on the left in both Figures 7 and 8, was chosen due to its low population level.The visual inspection carried out a few hours after the imagery acquisition confirmed the existence of a very small number of bees partially covering only one frame.The survival chances of colonies like this are very limited.In this case, all bees had disappeared a few days after the data acquisition.The second hive, shown in the middle image featured in Figures 7  and 8, was selected as a representative for those hives where routine monitoring is necessary.Their population level is low, partially covering less than a half of the brood chamber stacks.If their population level continues decreasing, maintenance tasks would be required to try to prevent the extinction of the colony.The third hive, shown on the right in Figures 7 and 8, was selected as a representation of strong colonies with very high populations.In this case, almost all frames are fully covered by bees producing large crowds at the top of the stacks due to the high density of bees inside the hive.In this case, the hive does not need immediate attention, being able to face the wintering phase successfully.Data acquisition was performed in November, immediately after honey gathering works and the preparation of hives for the start of the wintering phase.This preparation consists of the removal of all supplementary stacks keeping only the "brood chamber" (bottom stack).Each frame of the "brood chamber" was examined in order to ensure that all hives have some stacks with breeding and the remaining are full of honey, ensuring the existence of food for the survival of the hive until the next beekeeping season.
Data acquisition was carried out in the morning before sunrise.The ambient temperature at the time of data acquisition was 10 • C.
Data acquisition was planned independently for each sensor according to the directives described in Section 2.2.1.Considering that the objective of this step is to achieve a high-resolution 4D point cloud, the GSD was fixed to 1 mm for the RGB sensor and 5 mm for the thermographic sensor allowing the performance of image acquisition at a distance from the hive of up to 14 m and 2 m, respectively.
The RGB images were processed according to the photogrammetric and computer vision methodologies described in Section 2.2.2, obtaining a point cloud of each hive with an average resolution of more than 1,000,000 points/m 2 .The geometric accuracy of the reconstructed 3D RGB cloud was validated by comparing the main dimensions of the hives (i.e., width, height, and diagonals of each hive face) between the reconstructed 3D RGB cloud and the real hive.As a result, an averaged deviation of 3.5 mm with a maximum deviation of 5 mm was obtained.In addition, point-to-plane deviations were used to validate the reconstruction noise.Assuming that each hive face is a perfect plane, the least squares algorithm was used to compute and analyze point-to-plane deviations.These deviations present a normal distribution over the surface of the hive, with a mean of 0.8 mm and a standard deviation of 1.3 mm.The thermographic texture was projected over the point cloud obtaining a hybrid 4D point cloud with thermographic texture mapped over the points of the hive.This 4D point cloud with thermographic texture is subjected to a manual segmentation process to remove those parts that are not of interest, such as the top cover, floor, legs, or surrounding elements that do not belong to the hive, as shown in Figure 9. of 0.8 mm and a standard deviation of 1.3 mm.The thermographic texture was projected over the point cloud obtaining a hybrid 4D point cloud with thermographic texture mapped over the points of the hive.
This 4D point cloud with thermographic texture is subjected to a manual segmentation process to remove those parts that are not of interest, such as the top cover, floor, legs, or surrounding elements that do not belong to the hive, as shown in Figure 9.The inner temperature of the hive was calculated at each 3D point.To this end, information about the material and thickness of the hive walls was required, as shown in Equation (6).A digital superficial thermometer (Testo 720) was used to measure the interior temperature at certain points inside the hive.The comparison of the thermal values of the camera and the contact thermometer were inside the confidence interval defined by the contact thermometer manufacturer of ±0.2°.So, the confidence level of the results is limited by the accuracy of the contact thermometer.In this case, hive walls were built with wood with a thermal conductivity U of 0.13 m K/W and a thickness of 0.03 The inner temperature of the hive was calculated at each 3D point.To this end, information about the material and thickness of the hive walls was required, as shown in Equation (6).A digital superficial thermometer (Testo 720) was used to measure the interior temperature at certain points inside the hive.The comparison of the thermal values of the camera and the contact thermometer were inside the confidence interval defined by the contact thermometer manufacturer of ±0.2 • .So, the confidence level of the results is limited by the accuracy of the contact thermometer.In this case, hive walls were built with wood with a thermal conductivity U of 0.13 m K/W and a thickness of 0.03 m.Hive walls are topped with a thin layer of paint; however, the effect of the paint on the calculation can be dismissed due to its negligible thickness and to the fact that it would only affect the emissivity value of the surface.
The 4D model with thermal information of the interior of the hive is subjected to a pass-through filter to evaluate the surface in the hive interior that is over 12 • C, see Figure 10.A normal behavior of the temperature distribution on the inner surface of the hive can be seen.Following the upward flow of hot air, the surface temperature is higher in the upper area of the stack.A direct relationship between the inner surface of the hive walls that the colony is able to maintain above 12 • C and the population density can be seen.m.Hive walls are topped with a thin layer of paint; however, the effect of the paint on the calculation can be dismissed due to its negligible thickness and to the fact that it would only affect the emissivity value of the surface.
The 4D model with thermal information of the interior of the hive is subjected to a pass-through filter to evaluate the surface in the hive interior that is over 12 °C, see Figure 10.A normal behavior of the temperature distribution on the inner surface of the hive can be seen.Following the upward flow of hot air, the surface temperature is higher in the upper area of the stack.A direct relationship between the inner surface of the hive walls that the colony is able to maintain above 12 °C and the population density can be seen.

Conclusions
This article presents a non-invasive methodology for the analysis of the population condition and growing capabilities of honeybee hives through the generation of 4D models with geometric (X,Y,Z) and inner temperature (T) information.An automatic process for the analysis of the survival options of colonies for the wintering phase based on the evaluation of the colony capability to

Conclusions
This article presents a non-invasive methodology for the analysis of the population condition and growing capabilities of honeybee hives through the generation of 4D models with geometric (X,Y,Z) and inner temperature (T) information.An automatic process for the analysis of the survival options of colonies for the wintering phase based on the evaluation of the colony capability to maintain the inner wall temperature over 12 • C is proposed.
The described methodology can be used both in the research field for future evaluations of the behavior of hives and understanding of their survival processes and threats, as well as for the search of correlations towards the quantification of the population in the hive.Regarding the production field, the methodology can be used as a predictive tool to help in the decision-making of maintenance tasks.This methodology presents a more rigorous evaluation technique compared to other colony population assessments, with application regardless of the construction of the beehive and in a non-invasive way, since the opening of the beehive is not needed.
A deeper analysis will be performed in order to validate the methodology in hives made from different materials and formed of different shapes towards the knowledge of the distinct evolution of the colonies for the different hive configurations.In addition, long-term monitoring of the hives will improve the computation methodology of heat transfer, enabling the establishment of equations for the quantification of the population, and the identification of the most important heat transfer method for each time of the year.Authors are planning a complementary study where different beehives, placed in distant locations and under different conditions are analyzed through different techniques in different epochs.The study will compare the results and propose the optimum non-invasive beehive monitoring methodology, validating and probably improving the methodology presented in this work.
The generation of a high resolution and accurate 4D model not only allows for the selective measurement of temperatures on the points of interest but also the automation of analytical processes (segmentation based on critical temperatures).This work also opens new trends for future research allowing accurate monitoring of the hives through time to detect the behavior of the colonies against different pathologies that can cause their extinction.The characterization of the thermal behavior of the reduction in the population of the colony, either produced by external or internal agents, allows for scheduling of the necessary maintenance tasks that would avoid both the negative impacts in the internal conditions of the hive.In addition, the non-invasive character of the methodology would eliminate the stress suffered by the colony during the visual inspection inside the hive during the wintering phase.

Figure 1 .
Figure 1.Example of the images generated by the different imaging sensors used.(Left) Red Green Blue (RGB) image.(Right) Example of the thermographic image acquired, with thermal values represented using a color map: Dark blue is given to the lowest values (in this case, 2 °C), going to light blue and green for the highest temperature values (14 °C in this case study).

Figure 1 .
Figure 1.Example of the images generated by the different imaging sensors used.(Left) Red Green Blue (RGB) image.(Right) Example of the thermographic image acquired, with thermal values represented using a color map: Dark blue is given to the lowest values (in this case, 2 • C), going to light blue and green for the highest temperature values (14 • C in this case study).

Figure 2 .
Figure 2. "Master" and "slave" images following a "ring" in the data acquisition protocol.

Figure 2 .
Figure 2. "Master" and "slave" images following a "ring" in the data acquisition protocol.

Figure 3 .
Figure 3. 4D point cloud with thermographic texture of one beehive.Front and right sides (Left).Back and left sides (Right).

Figure 4 .
Figure 4. Components of incident InfraRed (IR) radiation acquired by the thermographic camera.

Figure 3 .
Figure 3. 4D point cloud with thermographic texture of one beehive.Front and right sides (Left).Back and left sides (Right).

Figure 3 .
Figure 3. 4D point cloud with thermographic texture of one beehive.Front and right sides (Left).Back and left sides (Right).

Figure 4 .
Figure 4. Components of incident InfraRed (IR) radiation acquired by the thermographic camera.

Figure 4 .
Figure 4. Components of incident InfraRed (IR) radiation acquired by the thermographic camera.

Figure 5 .
Figure 5. Diagram of the energy balance performed for the calculation of the inner temperature of the beehive: Heat flowing through the hive wall by conduction is equal to the heat exchanged from the outer surface of the hive wall to the outer air by convection plus the heat exchanged by this surface by radiation; both sources of radiation were received by the thermographic camera.

Figure 5 .
Figure 5. Diagram of the energy balance performed for the calculation of the inner temperature of the beehive: Heat flowing through the hive wall by conduction is equal to the heat exchanged from the outer surface of the hive wall to the outer air by convection plus the heat exchanged by this surface by radiation; both sources of radiation were received by the thermographic camera.

Figure 7 .
Figure 7. Visual inspection from a top view to determine the population through the count of frames full of bees (highlighted in red).Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 7 .
Figure 7. Visual inspection from a top view to determine the population through the count of frames full of bees (highlighted in red).Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 7 .
Figure 7. Visual inspection from a top view to determine the population through the count of frames full of bees (highlighted in red).Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 8 .
Figure 8. 4D models which integrate metric and thermographic information.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 8 .
Figure 8. 4D models which integrate metric and thermographic information.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 9 .
Figure 9. Segmentation of the hives: Surfaces to be evaluated of the three hives for their population assessment.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 9 .
Figure 9. Segmentation of the hives: Surfaces to be evaluated of the three hives for their population assessment.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 10 .
Figure 10.Surface evaluation over pass-through filter results.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Figure 10 .
Figure 10.Surface evaluation over pass-through filter results.Very low level of population (Left).Mid level of population (Middle).High level of population (Right).

Table 1 .
Interior orientation parameters of visible camera Canon EOS 700D equipped with a 60 mm fixed focal lens, as a result of its geometric calibration.

Table 1 .
Interior orientation parameters of visible camera Canon EOS 700D equipped with a 60 mm fixed focal lens, as a result of its geometric calibration.

Table 2 .
Interior orientation parameters of thermographic camera, NEC TH9260, resulting from its geometric calibration.