Mapping almond orchard canopy volume, flowers, fruit and yield using lidar and vision sensors

14 Pages·2016·6.73 MB·English

by  James P. Underwood

ComputersandElectronicsinAgriculture130(2016)83–96
ContentslistsavailableatScienceDirect
Computers and Electronics in Agriculture
journal homepage: www.elsevier.com/locate/compag
Original papers
Mapping almond orchard canopy volume, flowers, fruit and yield using
lidar and vision sensors
⇑
James P. Underwood , Calvin Hung, Brett Whelan, Salah Sukkarieh
AustralianCentreforFieldRobotics,TheUniversityofSydney,NSW2006,Australia
a r t i c l e i n f o a b s t r a c t
Articlehistory: Thispaperpresentamobileterrestrialscanningsystemforalmondorchards,thatisabletoefficiently
Received15April2016 mapflowerandfruitdistributionsandtoestimateandpredictyieldforindividualtrees.Amobilerobotic
Receivedinrevisedform20September groundvehiclescanstheorchardwhileloggingdatafromon-boardlidarandcamerasensors.Anauto-
2016
matedsoftwarepipelineprocessesthedataoffline,toproducea3Dmapoftheorchardandtoautomat-
Accepted24September2016
icallydetecteachtreewithinthatmap,includingcorrectassociationsforthesametreesseenonprior
occasions.Colourimagesarealsoassociatedtoeachtree,leadingtoadatabaseofimagesandcanopy
models,atdifferenttimesthroughouttheseasonandspanningmultipleyears.Acanopyvolumemeasure
Keywords:
isderivedfromthe3Dmodels,andclassificationisperformedontheimagestoestimateflowerandfruit
Robotics
density.Thesemeasureswerecomparedtoindividualtreeharvestweightstoassesstherelationshipto
Sensing
Machinevision yield.Ablockofapproximately580treeswasscannedatpeakbloom,fruit-setandjustbeforeharvestfor
Lidar twosubsequentyears,withupto50treesindividuallyharvestedforcomparison.Lidarcanopyvolume
Multi-sensorfusion hadthestrongest linear relationshipto yield withR2¼0:77for 39 treesamples spanning twoyears.
Orchardyieldmapping Anadditionalexperimentwasperformedusinghand-heldphotographyandimageprocessingtomeasure
fruitdensity,whichexhibitedsimilarperformance(R2¼0:71).Flowerdensitymeasurementswerenot
stronglyrelatedtoyield,however,themapsshowcleardifferentiationbetweenalmondvarietiesand
maybeusefulforotherstudies.
(cid:1)2016ElsevierB.V.Allrightsreserved.
1.Introduction Thesystemcontinuouslyrecordsdatafromacameraandlidarsen-
sor,whilethevehicledrivesthroughtheorchard.Thisiscombined
Technologicalimprovementsinsensing,computing,algorithms with processing software that automatically extracts geometric
and robotics have thepotential toincrease productivityfor com- and visual information of each tree and matches the data from
mercial farming and efficiency for plant science. For farmers, scans taken at different times of the year and over multiple sea-
mobile data gathering systems can provide detailed information sons.Thisallowstheassessmentoffloweringandfruitproduction
to assist their decision making and management processes and perindividualtree,inamannerthatisefficientforscanningwhole
the information can plug into decision support software that is orchardblocks.Thisrepresentsasignificantincreaseinresolution
capableofrecommendingparticularactions.Eventuallyitwillbe comparedtothetypicalpracticeofweighingthetotalproducehar-
possiblefortheseactionstobedirectlyappliedusingmobilefield vested from whole orchard blocks, allowing the variability
roboticstechnology.Forplantscientists,mobiledatasystemscan betweenindividualtreesandsmallerregionsoftheorchardtobe
providehighthroughput,in-fieldplantphenomics.Thiswillallow estimatedandmapped.
greatercapacityforin-fieldexperimentation,wheremanuallabour Threedimensionalandimage-basedsensinghavebeenapplied
forin-fielddataacquisitioniscurrentlyalimitingfactor,leadingto tomanyaspectsoftree-cropprecisionagriculture.Therearemany
yield improvements from genomics and improvements to best- examples of the use of lidar to measure tree canopy geometry
practiceforgrowers. (Tumbo et al., 2002; Walklate et al., 2002; Rosell et al., 2009;
This paper presents a robotic ground-vehicle information sys- Rosell and Sanzs, 2012; Wellington and Campoy, 2012; Escolà
tem for almond orchard mapping and per-tree yield estimation. et al., 2015) and as a proxy for related measurements such as
leaf-area-index(Sanzetal.,2013).Alternativerangesensorshave
⇑ also been used including ultrasound (Tumbo et al., 2002;
Correspondingauthor.
Hosainpouretal.,2013),structuredlight(Rosell-Poloetal.,2015)
E-mail addresses: [email protected] (J.P. Underwood), brett.
andstereovisionRosellandSanzs(2012),butlidarispopulargiven
[email protected](B.Whelan),[email protected](S.Sukkarieh).
http://dx.doi.org/10.1016/j.compag.2016.09.014
0168-1699/(cid:1)2016ElsevierB.V.Allrightsreserved.
84 J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96
therelativelyhighaccuracyandinvarianceundernaturalillumina- grape berry detection method in Nuske et al. (2014). Algorithms
tionconditions. are also tailored to the orchard configuration, such as the two
Vision sensors have been coupled with machine vision algo- dimensional trellis fruit wall in Hung et al. (2015), which avoids
rithms to estimate fruit and flower densities for individual trees the need for individual tree segmentation for image masking,
for a number of different crop types (Gongal et al., 2015). Hand- whichisakeycomponentofoursystem.
held digital cameras and relatively simple image classification Although several of the components that are necessary for an
algorithms have been used to estimate flower densities almondorchardscanningsystemhavebeenexploredinthelitera-
(Adamsen et al., 2000; Aggelopoulou et al., 2011; Thorp and ture(includingtheauthors’priorworkandothers),thecontribu-
Dierig, 2011); relatively simple algorithms are possible due to tion of this paper is to develop an integrated methodology to
the typically high colour contrast exhibited by flowers. Machine create a single, efficient orchard scanning system for almonds.
vision cameras have also been mounted on tractors to improve Thecontributionisthecompletesystem,includingthenecessary
theimageacquisitionprocess(Hocˇevaretal.,2014),whichallowed developments to combine the components of tree detection and
flowerestimationonlargernumbersoftrees(N=136).Neverthe- segmentationwithflowerandfruitmappingandyieldestimation,
less, the process requires manual demarcation of the individual together with experimentation covering approximately 580 trees
trees withinthe frames,whichlimitsthe efficiencywhenscaling atthreetimesofseasonfortwoyears.
uptoscanentirecommercialorchardblocks.Unlikeflowers,fruit
andnutsare oftenvisuallysimilar totheleavesandsurrounding
2.Materialsandmethods
foliage of trees, meaning that more sophisticated machine vision
methodsarerequiredtoautomaticallydetectthem.Theseinclude A 2:3-hectare block of a commercial almond orchard was
artificial neural networks (Gongal et al., 2015; Hung et al., 2013,
scanned using the ‘‘Shrimp” ground vehicle robot, three times
2015;BargotiandUnderwood,2016a,b),whichhavetheadvantage
per year (flowering, fruit-set, pre-harvest) for two subsequent
ofautomaticallylearningappropriatefeaturedescriptionsforclas-
years.Atharvest,individualtreeswerealsophotographedwitha
sificationfromthedata;multi-sensorfusion,whichsimplifiesthe
hand-held digital camera and then selectively harvested and
classificationproblembyfusingdatafromcomplementarysensors
weighed.Thecolourimagesandlidarfromtherobotandtheman-
suchasvisionandthermalcameras(Bulanonetal.,2009)orvision
uallytakenphotographswerepost-processed,usingcustomalgo-
and near-infra-red cameras (Hung et al., 2013); and hybrid
rithms and software, to derive measures relating to the canopy
approaches that combine colour and shape (Singh et al., 2010;
volume and the density of flowers and fruit on each tree. The
Nuske et al., 2011; Wang et al., 2013). Cameras have also been
measures were compared with the selective harvest weights, to
combined with lidar for tree-crop applications, such as Shalal
quantify performance. This section describes the system and
et al. (2015) and Bargoti et al. (2015), which both address tree
sensors, the protocol for data collection, and how the data were
trunk detection by combining the geometry sensed by the laser
processed.
withthevisualappearancesensedbythecamera.
Thesizeofalmondtreesisknowntobeanimportantfactorin
estimating the yield (Hill et al., 1987), which motivates canopy 2.1.The‘‘Shrimp”mobilerobot
geometrysensing.Flowerdensitiesarealsoconsideredtoberele-
vanttoyield,althoughtherelationshipiscomplicatedbyvariabil- The ‘‘Shrimp” mobile ground vehicle robot was designed and
ity in pollination (e.g. availability of pollinators) and other built at the Australian Centre for Field Robotics at the University
limitations in how much fruit the tree can bare Dicenta et al. ofSydneyin2009(seeFig.1(a))asageneralpurposeresearchplat-
(2005). The potential utility of flower density mapping, as well formtostudyroboticsensingandperception.Forthisapplication,
potentiallybeingabletodirectlymeasurefruitdensitymotivates weuseasubsetofsensors:a2Dlinescanninglidar(SICKLMS-291)
theuseofvision. and a machine vision camera (Point Grey Ladybug2, single two
Inordertoallowallofthesemethodstoscaleuptoentireorch- mega-pixel camera, with natural illumination only), both face to
ards, automated, streamlined data management is also required, therighttoscanthetreesasthevehicletravelscontinuouslyfor-
which includes software for tree segmentation and detection wardsat upto 2m/s(see Fig. 1(b)). Areal timekinematicglobal
(Wellingtonand Campoy,2012;Shalalet al.,2015;Bargotietal., positioning inertial navigation system (Novatel SPAN RTK
2015; Underwood et al., 2015b) and tree matching (correct data GPS/INS) is used for positioning, a gamma radiometer (RS700)
to tree assignment) for repeated scans at different times was mounted on-board to record passive soil gamma emissions
(Underwoodetal.,2015b).Therearefewwhole-systemexamples and an electromagnetic induction instrument (EM38) is towed
thatcombinegeometricandvisualsensing,togetherwithefficient behindthevehicletomeasureapparentsoilelectricalconductivity
mobiledataacquisitionandautomateddataprocessingandman- (ECa).Thesystemincludesacomputerfordatalogging.
agement steps that facilitate entire blocks of commercial farms
to be efficiently scanned, including comparisons to ground truth 2.2.Datacollection
yieldsuchasfruitcountsorharvestweights.Prominentexamples
includeholisticsystemsforgroundcrops(Busemeyeretal.,2013), Alldatawereobtainedfroma2.3-hectaresectionofLakeCul-
vineyards (Nuske et al., 2014) and apple orchards (Hung et al., lulleraineAlmonds,inVictoria,Australia,showninFig.2.Thearea
2015). Each of these systems relies on constraints relating to the includes 10 rows spaced 7.35m apart and roughly 313.5m long,
specificnatureofthetargetcropandnoonesystemandapproach with58treesperrowspacedat5.5m.Theprimaryvarietyinevery
islikelytobeadaptabletovastlydifferentcroptypes.Thegeome- second row is Nonpareil, with alternating Carmel and Monterey
try of ground based crops is well suited to systems that straddle pollinator rows between. Datasets were collected at peak bloom
above the crop (Busemeyer et al., 2013; Deery et al., 2014; (asestimatedbythefarmmanager),atfruit-setandjustpriorto
Underwood et al., 2015a), which enables controlled illumination harvest,forboththe2014and2015harvestseasons.Severaldata-
andprovidesforanidealsensorvantagepoint,butthosesystems sets were taken on subsequent days to assess repeatability. All
are not applicable to tree crops, which are taller and difficult to datasetsaresummarisedinTable1.
straddle.Amongsttreecropapplications,algorithmsaretypically To obtain each dataset, the raw data from all sensors were
tailoredtotheappearanceandgeometryofthespecificfruitsuch loggedcontinuously,whilethevehicleisdrivenbyaremoteoper-
ascircledetectionforapplesinHungetal.(2015)andacustomised atorupandbackdowntherowsatspeedsof1–2m/sdependingon
J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96 85
Fig.1. (a)The‘‘Shrimp”mobilerobot,equippedwithavarietyofsensors,includingacolourcamera,side-scanninglidar,EM38soilconductivity,RS700gammaradiometer
andGPS/INSlocalisationsystem.(b)Theverticalcurtainoflidardataissweptacrossthetreesasthevehiclemovesforward(totheleftoftheimage).(c)Atypicalimageofa
singlealmondtreeatfruit-set.
Fig.2. LakeCullulleraineAlmonds,Victoria,Australia(a)andlabelledimage(b)showingthe2.31hablockthatwasscannedinthisstudy.Thepinin(a)denotesthetrailer
basestationmarkedin(b).Rowsareseparatedby7.25mandtreetrunksbyanaverageof5.5m.SatellitephotosfromGoogle2013.
Table1
conditions. The sensors face to one side of the vehicle, so each
Datasets.
access row is traversed in both directions to scan both sides of
Year Date Season thetrees.Itisimportanttoscanbothsidesbecauseocclusionspre-
1 2013-08-15 Flowering(day1) ventthewholetreefrombeingcaptured(withlidarorvision)from
2013-08-16 Flowering(day2) onlyoneview-point(LinandHyyppä,2012).Inthefirstyear,the
2013-09-23 Fruit-set
vehicle was driven up the first row, performed a u-turn at the
2014-02-10 Pre-harvest
end and then back down the same row (sensors facing the other
2 2014-08-26 Flowering(day1)
way), repeated for all rows in sequence. In the second year, the
2014-08-27 Flowering(day2)
2014-10-14 Fruit-set(day1) vehiclewasdrivenforwardsandbackwardsinanalternatingpat-
2014-10-15 Fruit-set(day2) tern,halfinthemorningfacingwestandhalfintheafternoonfac-
2015-01-28 Pre-harvest ingeast,toensurethesunwasalwaysbehindthecamera,toavoid
86 J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96
Row of Trees
Fig.3. (a)Theprotocolforphotographingindividualtreeswithahand-heldcamera.Thelargeandsmallcirclesrepresentcanopiesandtrunksrespectively.Thetargettreeis
showningrey,fourphotosaretakenfromthebaseofthearrowspointingtowardsthetrunkofthetargettree.(b)Atypicalphoto.
lens-flares and saturation that affected the data in the first year. 2.3.Dataprocessing
Thetotaltraversalforeachscanwas6.2km,takinganaverageof
1.5h.1 The lidar data was calibrated and georeferenced with the
Atharvesttime,atotalof30individualtrees(11Nonpareil,9 GPS/INS localisation system (Underwood et al., 2010) to produce
Monterey,10 Carmel) in thefirst year and 50 in the second year a 3D point cloud of the entire scanned area. Individual rows are
(20 additional Nonpareil) were selectively harvested. The total automatically separated within the data by detecting where the
numberwasconstrainedbytime,labourcostandminimisingdis- vehicleturnsattheendsoftherow.Individualtreesaresegmented
ruptiontothecommercialorchard.Theinitialsamplingpatternof using the technique first established by Wellington and Campoy
30 trees was selected by balanced random sampling within five (2012)andextendedbyUnderwoodetal.(2015b),asdepictedin
classes of observed soil variation. The classes were created using Fig.1(b),wheredatafromeachtreeiscoloureddifferentlytoindi-
a multivariate hierarchical clustering process following the catethesegmentation.Theresultisthatdatarelatingtoeachtree
method of Taylor et al. (2007) with soil ECa and the total count canbeuniquelyidentifiedandseparated,accordingtoitsrowand
of soil gamma emissions obtained in the initial survey used as tree number and compared between datasets taken at different
inputlayers.Thefieldwasstratifiedusingthesensedsoildatato times.Thelidarpointcloudofeachtreeisusedtocalculateamea-
ensureanyproductionvariationdrivenbyobservedvariabilityin sureofthecanopyvolume,describedinSection2.3.1.Asinglecor-
soilpropertieswascapturedintheyieldsampling.Thepatternof respondingphotocentredoneachtreecanbeextractedfromthe
20additionaltreeswaschosentoexpandtheNonpareilcountwith freerunningvideodataautomatically,3asillustratedinFig.4.The
auniformspreadoflidarcanopyvolume,giventherelationshipto georeferenced lidar data describe the canopy geometry precisely.
yieldseenfromthe2013data(describedinSection3).Theground Thedataforeachindividualtreeareprojectedintothecorrespond-
underneath each tree was swept clear of debris, the trees were ing image via a process of camera-lidar-platform calibration and
mechanically shaken and the nuts were collected and weighed. geometric transformation (Underwood et al., 2007). The result is
Inthefirstyear,acrack-outwasperformedbytaking1–2kgsam- an image overlay that marks the position of the canopy as shown
ples,crackingandweighing100kernelstoestimatethetotalnum- inFig.4(d),whichisusedasaregionofinterestforfurtherprocess-
berofkernels,averagekernelweightandtotalkernelweightper ing to estimate flower and fruit abundance, as described in Sec-
tree. In the second year, only the total harvest weights per tree tion 2.3.2. Measures relating to canopy volume, flower and fruit
wereabletobeobtained. densityarederivedfortwofullseasons.Thesearemappedandcom-
Inthesecondyear,ahand-heldcamera(CanonEOS-60D,17.9 pared to the post harvest weights for the individually harvested
mega-pixel2) was used to manually photograph the 50 sampled trees,toquantifytheperformanceofthesystem.4Theentireprocess
trees,tocomparetheperformanceofalmonddetectionwithahigher was run offline from logged sensor data. Precise processing times
resolutionimage,manualcontroloflightingconditionsandmultiple were not recorded, however, the lidar processing component took
perspectivespertree.Asimpleprotocolwasusedtotakeconsistent less than 5min on a typical computer, to segmentall 580 trees in
photographs, illustrated in Fig. 3. Four photos were taken of each the2.3-hectareorchardblock.Themosttimeconsumingcomponent
tree, from diagonal positions to maximise the distance from the (includingbothprocessingtimeandhumanlabellingeffort)wasfor
camera to the tree, while avoiding occlusion from neighbouring theimageprocessingstepsdescribedinSection2.3.2,whichrequire
trees.Allphotosaretakenwiththesunbehindthecamera,achieved approximatelytwodaysofcombinedhumanandcomputertimeper
by taking west/east facing photos in the morning/evening 2.3-hectaredataset.
respectively.
2.3.1.Treecanopyvolumefromlidar
The3Dpointcloudmodelofeachtreeisprocessedseparately.
1 Ifduplicatesensorsfacedleftandright,thetrajectoryandscantimecouldbe Althoughthecanopiesareapproximatelyellipsoidal,theyexhibit
halved,butcarewouldberequiredtoscanasclosetothezenithaspossible.Both
trajectoriesdescribedherehavethesametotallength.
2 Thehighestresolutioncamerathatwasavailabletousfortheseexperimentswas 3 Fromthelidar,weknowexactlywhereeachtreeisandwhenitwasscannedby
chosen,withanimagesizeof3456by5184pixels,tomaximisethepixel-countper Shrimp,sowecanextractcorrespondingdatafromallothersensorsautomatically.
almond. 4 Avideoillustratingtheprocessisavailablehere:http://youtu.be/j-gJQXQoqX0.
J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96 87
Fig.4. Photographsofthesamealmondtreeduring(a)flowering,(b)fruit-setand(c)justpriortoharvest.Thephotoswereautomaticallysegmentedfromvideotakenbythe
robot‘‘Shrimp”,byco-registrationtothelidarmap.Thelidarcanopygeometryisprojectedintotheundistortedimage(d)asaregionofinterestforfurtherprocessing.
Fig.5. A3Dpointcloudmodelofanindividualtreeaftersegmentation.Thetreeisscannedfrombothsidesatdifferenttimes,shownbydifferentcolours.Awell-aligned
modelin(a)perspectiveviewand(b)sideviewand(c)amisalignedmodelduetolocalisation(GPS)error.Misalignmentcansignificantlyaffectthecalculationofcanopy
volume.(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
amorecomplexinternalstructure,withbranchesandvoidspace For this tree, the lidar estimates a foliage volume of about 45m3,
within. Much of this structure is detectable by the lidar scans, whichisexpectedgiventhatalargeamountoftheinternalcanopy
becausethelaserspartiallypenetratethroughtheouterlayersof volume is free space. Each individual tree is automatically seg-
the tree. We therefore model the canopy foliage volume by dis- mentedbyoursoftwarepipeline,soweareabletoefficientlycalcu-
cretisingthe points intocubic voxels and tallying the total num- late the canopy volume per tree. Our method avoids the need to
ber.5 Multiplying the count by the voxel size (0.001m3) gives an align the data, however, the problem of aligning lidar scans from
estimateofthecanopyfoliagevolume.6Becauseeachsideofatree bothsidesofthetreehasbeenstudiedintheliterature,usingman-
is scanned separately, errors in the platform localisation that are ualalignment(Roselletal.,2009),whichmightnotefficientlyscale
caused by GPS drift mean that the two tree halves cannot always tohundredsoftreesandusingsimultaneouslocalisationandmap-
be accurately aligned using the localisation data, as illustrated in ping(SLAM)(CheeinandGuivant,2014),whichshouldscalewell.7
Fig. 5, with an example of accurate alignment in Fig. 5(a) and (b),
andanexampleofmisalignmentinFig.5(c).Thesensitivityoftree
volume errors due to trajectory was experimentally studied in 2.3.2.Flowerandfruitdensityfromimagery
Pallejaetal.(2010)andacomprehensivemathematicalandexperi- Flowerandfruitdensityestimatesareproducedbyperforming
mentalanalysisforlidaringeneralispresentedinUnderwoodetal. imageclassificationontheimagesassociatedtoeachtree,within
(2010,2009).Wethereforecalculatethevolumesofthetwohalves thelidarcanopymask.Shrimpwasusedtoscantheorchardblock
separately and sum them together. This eliminates the effect of duringflowering,fruit-setandpre-harvestperiodsfortwowhole
misalignment, but due to double counting the overlapped region, seasons,resultinginadatabaseofphotographsforthesekeytimes.
this may cause a bias depending on the geometry of the canopy. Fig.4showsthephotosofthesametreeretrievedfromthedata-
The well aligned canopy in Fig. 6 is 5.58m by 6.35m wide and base at flowering, fruit-set and pre-harvest in the first year. The
4.52mtall,forasolidellipsoidalcanopyvolumeofroughly84m3. 3D lidar mask (see Fig. 4(d)) indicates the location in the image
where the tree must be according to the laser scan, then fruit
andflowerdetectionareperformedwithinthatregionoftheimage
5 Aresolutionof10cmwasused,thoughtheoutputwasfoundtobeinsensitiveto
otherresolutions. 7 SLAMforcanopyvolumeestimation(CheeinandGuivant,2014)wouldtheoret-
6 Note that lidar returns from the ground are excluded from calculations by icallybeappropriatefortreecanopyalignmentinthecontextofoursystem,which
removingpointswithin20cmelevationofwherethewheelsofthevehiclepassed, wouldalsoallowanybiasesinourmethodtobequantified,butthiswasoutsidethe
however,thegroundisshowninthefiguresforcontextualisation. scopeofthisstudy.
88 J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96
Fig.6. A3Dvoxelisedpointcloudofatree,colouredbyelevation.Thecomplexstructurecanbeseenandispreservedinthecalculationofcanopyvolume.(Forinterpretation
ofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
to avoid confusion due to adjacent trees. For the hand-held pho- is more efficiently reduced by simple morphological operations
tographs, the lidar mask is unavailable, so the canopy masks are (erosionanddilation)insteadoftheCRFstep.
manuallydrawn. Anexampleofarawimageatfruit-setandtheresultingclassi-
Forbothflowerandfruitestimation,imageclassificationisper- ficationisshowninFig.7(b).Theclassifierisabletodiscernmuch
formed and the measured variable is the total number of pixels of visually larger fruit in the foreground, however, the algorithm
thatareclassifiedasfruitorflowerswithinthelidar-basedimage failstodetectsmallerfruitinthebackground,particularlywhere
mask, divided by the total number of pixels in the mask. This is individual almonds are silhouetted against the bright sky. Hand-
referredtoasflowerorfruitpixeldensity.Thisformulationisused held photos were also captured for 50 individual trees with a
to normalise the measurement with respect to variation in per- higherresolution,withanexampleshowninFig.8.Qualitatively,
spective.Thisisnecessarybecauseslightdifferencesintheproxim- it can be seen that the classifier is better able to detect fruit in
ityofthecameratothetreecanaltertheapparentsizeofthetree thehand-heldimages,althoughthereissomeconfusionbetween
intheimage.Thedensityestimateforbothflowersandfruitisa thewoodyalmondsjustbeforeharvestandthelimbsofthetree.
dimensionlessratioandnotadirectestimateoftheabsolutenum-
berofflowersorfruit.Evenifitwerepossibletocounttheabsolute
3.Results
numberofvisiblefruitorflowerswithintheimage,visualocclu-
sions mean that this quantity is only a fraction of the unknown
This section quantifies the performance of the different mea-
truecountonthewholetree.Wethereforeassumethattheratio
sures, including lidar canopy volume, flower and fruit density.
ofvisiblefruitorflowerstothetotalisconstantbetweendifferent
Where possible, the consistency of the measures is quantified in
trees.Specifically,weassumethatifNflowersorfruitareseenin
termsofrepeatability,whichisanecessary(thoughnotsufficient)
the imagesfrom the two opposing sides of a tree, then there are
conditionfor aninstrumentto provideaccurateinformation.The
N(cid:2)k flowers or fruit in total for some unknown constant k that
performance of all measures were qualitatively assessed against
is determined from calibration with actual yield data. Wherever
observations from thefield and quantifiedagainstthe selectively
this assumption does not hold well, we expect a non-linear or
harvested fruit weights, to assess the ability of the measures to
noisy relationship between the vision estimates and the ground
predictyield.
truth. This assumption is tested in this work by measuring the
repeatability and consistency of the relationship of image esti-
matestoactualyield. 3.1.Lidarvolumeconsistency
In our dataset the flowers could be separated easily from the
backgrounddue tothe highlycontrasting whitecolour,therefore The consistency of the lidar volume measure was assessed by
a simple colour threshold was sufficient: flower pixels were comparing the measurements of the same trees (approx. 580 in
defined as hue>0:51;saturation<0:3;value>0:9, with all total),usingdatasetsobtainedatdifferenttimeswithineachsea-
son.Inthefirstyear,thisincludedtwodatasetsduringflowering
parametersdeterminedempiricallyfromthedata.Theupperlimit
(onsubsequentdays),oneduringfruit-setandonejustbeforehar-
wasplacedonsaturationtoavoidfalselydetectingovercastskies
vest.Inthesecondyear,wecomparedtwosubsequentdaysduring
as flowers. An example image and flower detection output is
flowering, two subsequent days at fruit-set and one just before
showninFig.7(a),whereitcanbeseenthatmostoftheflowers
harvest. Scans taken one day apart can be used to estimate the
are detected, with some false detection along the sunlit edge of
repeatabilityofthesystem,whereasforscanscomparedatdiffer-
oneofthebranches.Thecolourthresholdingapproachwasfound
enttimesofyear,thetreegeometrychanges,butwestillexpecta
to be sufficient for our data, but a more sophisticated method
simplelinearrelationshiptohold.
mayberequiredtohandleilluminationvariationinthefuture.
TherelationshipbetweenmeasurementsisquantifiedinTables
Forfruitdetection,weusethemulti-scalefeaturelearningalgo-
rithmfromHungetal.(2013,2015),withmodificationstoimprove 2and3bythecoefficientofdetermination(R2).Strongagreement
the efficiency of the algorithm. The original implementation is seen for all scans separated by one day (R2¼0:96;0:99;0:94).
appliedConditionalRandomFields(CRF)tosmooththeoutputof The repeatability is quantified by focussing on all 580 trees
the multi-class labels, however, this paper focuses on binary scannedonedayapartduringfloweringinthefirstseason,which
fruit/non-fruitclassificationandthenoiseofthebinaryclasslabels is shown in Fig. 9(a). The corresponding coefficient of variation
J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96 89
Fig.7. Rawimagesandthecorrespondingflower(a)andfruit(b)detections.
Fig.8. Aphotographofanalmondtree,takenwithahand-heldcameraandtheresultofmanuallymaskingthetreeandautomaticallyclassifyingthefruit(whitepixels).
90 J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96
Table2
Lidarvolumeagreement(R2)forrepeatedscansduringfirstyear((cid:3)580samples).
Flowering(day1) Flowering(day2) Fruit-set Pre-harvest Dataset Date
1 0.96 0.90 0.73 Flowering(day1) 2013-08-15
1 0.89 0.72 Flowering(day2) 2013-08-16
1 0.82 Fruit-set 2013-09-23
1 Pre-harvest 2014-02-10
Table3
Lidarvolumeagreement(R2)forrepeatedscansduringsecondyear((cid:3)580samples).
Flower(day1) Flower(day2) Fruit-set(day1) Fruit-set(day2) Pre-harvest Dataset Date
1 0.99 0.90 0.89 0.44 Flowering(day1) 2014-08-26
1 0.89 0.89 0.44 Flowering(day2) 2014-08-27
1 0.94 0.51 Fruit-set(day1) 2014-10-14
1 0.49 Fruit-set(day2) 2014-10-15
1 Pre-harvest 2015-01-28
Canopy Foliage Volume (m3) R=0.98 Lidar Canopy Volumes: Year 1 vs Year 2
70 70
3m)
3 (
60 1 60
0
2
3 g
1 n
20 50 eri 50
08 / Flow
5 / 40 ng 40
Flowering: 1 2300 mcnaoornmnpteealrreeyil py Volume Duri 2300 MCNaoornnmpteaerrleeyil
o
n
a
10 C
ar 10
d
Li
0
0 10 20 30 40 50 60 70 80 00 10 20 30 40 50 60 70 80
Flowering: 16 / 08 / 2013 Lidar Canopy Volume During Flowering 2014 (m3)
(a) Volume one day apart (b) Volume one year apart
Fig.9. Lidarcanopyvolumemeasurementsbetweensuccessive(a)daysand(b)years.Pointsonthelineindicatenochange,pointstotherightindicategrowth.
(CV)is2:4%,whichfallswithintherangeoflessthanthreepercent The lidar canopy volumes are shown for both years (during
reportedbyWeiandSalyani(2005)forlaser foliagedensityesti- flowering)inFig.9(b).Measurementsonthelineofequalityindi-
mationforcitrustrees.8Thevariationislikelyduetominordiffer- cate a tree with the same measured volume in both years. Most
ences in the vehicle trajectories (velocity, distance and viewing measurementsfalltotherightoftheline,indicatinggrowth.The
angle to tree)perturbingthe samplingpatternin the point clouds. smaller,replantedtreesincreasedinvolumemorethanthemature
Nevertheless, the canopy volume measurements are highly trees.CanopygrowthpervarietyisemphasisedinFig.10.Thedata
repeatable. reveal that on average, all trees increased their canopy size
Decreasing measurement agreement was observed with betweensuccessiveyears,buttheNonpareiltreesincreasedtheir
increasing temporal separation between scans, indicating that volumebymorethanthetwopollinatorvarieties.
the changes in canopy geometry are not perfectly predictable. E.
g.thetreevolumedoesnotuniformlychangebyafixedpercentage
3.2.Lidarvolumeandyield
ofthevolumemeasuredearlierintheseason.Thepoorestagree-
mentwasobservedwiththemeasurementstakenjustbeforehar-
Thelidarcanopyfoliagevolumeswerecomparedtothethe30
vest (R2 as low as 0.44). This was due to a degradation in the
individualtreesamplesfromthefirstyear(11Nonpareil,10Car-
qualityofrawdatawhenthefoliageismostdense:theGPSsatel-
meland9Monterey)includingkernelcountsandweightsfroma
litesignalsbecomemoreheavilyoccluded,thelidarislessableto
crack-outprocess.Inthesecondyear,comparisonsweremadeto
penetratethecanopy,andtheneighbouringcanopiesoverlapmore
atotalof31Nonpareiltrees(includingthesame11fromthefirst
heavily,whichchallengesthetreesegmentationprocess.
year),butonlytotalharvestweightwasabletobemeasured.Alin-
earregressionwasperformedandthecoefficientofdetermination
is reported. The results for the first year are shown in Table 4. A
8 Duetotheefficiencyaffordedbyoursystemwehaveincreasedthenumberof
highdegreeofagreement(R2>0:70)isobservedbetweenthevol-
samplesfrom10treesto580forourrepeatabilityestimates,comparedtoWeiand
Salyani(2005). umes and yield for all varieties, though it should be emphasised
J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96 91
Histogram of Canopy Growth Year 1 to Year 2 Histogram of Canopy Growth Year 1 to Year 2 Histogram of Canopy Growth Year 1 to Year 2
35 9 14
8
30 12
7
Number of trees122505 Number of trees 456 Number of trees1680
3
10 4
2
5 2
1
0 0 0
−10 −5 0 5 10 15 20 25 −10 −5 0 5 10 15 20 −10 −5 0 5 10 15
Growth (m3) Growth (m3) Growth (m3)
Fig.10. Histogramoflidarcanopyvolumegrowthpertreefromthefirsttosecondyearforthethreevarieties.Meangrowthislistedundereachfigure.
Table4
Canopyfoliagevolumerelationshipwithyield(R2,10samplespervariety).
Variety Totalweight Kernelweight Kernelcount Dataset Date
Nonpareil 0.83 0.70 0.65 Flowering(day2) 2013-08-16
0.77 0.58 0.52 Fruit-set 2013-09-23
0.64 0.38 0.31 Pre-harvest 2014-02-10
Monterey 0.77 0.80 0.74 Flowering(day2) 2013-08-16
0.73 0.78 0.75 Fruit-set 2013-09-23
0.50 0.58 0.56 Pre-harvest 2014-02-10
Carmel 0.88 0.80 0.82 Flowering(day2) 2013-08-16
0.66 0.60 0.61 Fruit-set 2013-09-23
0.61 0.58 0.64 Pre-harvest 2014-02-10
Averagea 0.82 0.77 0.74 Flowering(day2) 2013–08-16
0.72 0.66 0.63 Fruit-set 2013-09-23
0.58 0.51 0.50 Pre-harvest 2014-02-10
Combinedb 0.28 0.35 0.39 Flowering(day2) 2013-08-16
0.19 0.25 0.29 Fruit-set 2013-09-23
0.13 0.17 0.21 Pre-harvest 2014-02-10
a AverageR2overvarieties.Thisisnotastandardstatisticalterm,butisusedheretosummarisetheperformanceinthistable.Thesizeofthedatasetshouldbeconsidered
asN=10,notN=30becausethevarietiesaretreatedseparately.
b CombinedassumesallvarietiescanbelumpedtogetherinonebatchwithN=30.Thisshowsthatalthoughthereisastrongrelationshipbetweencanopyfoliagevolume
andyield,thelinearrelationshipisdifferentforeachvarietyandmustbemodelledseparately.
thatN(cid:3)10treesforeachvarietyissmall.9Ifallthreevarietiesare increasedtheiryieldbyapproximately175%and2000%,andthat
lumped together with a single linear relationship from volume to the most productive sampled tree produced approximately 150%
yield,theagreementispoor(R2<0:39).Thedatafromthefirstyear more yield than any other, which are all unlikely occurrences. It
suggeststhatcanopyfoliagevolumeagreeswellwithyieldbutthat would be extremely unlikely that our system would identify all
thelinearrelationship isdifferent foreachvariety.For thisreason, three cases (and no others) as outliers by random chance. The
samplinginthesecondyearfocussedonthemainNonpareilvariety, selective harvest was performed by low skilled farm labourers
toincreasethesignificanceoftheresults. whoaretrainedonlytoharvestasefficientlyaspossibleinnormal
In the second year, with 31 Nonpareil samples we observed a productionconditionsandnottometiculouslygatherdataforsci-
coefficientofdeterminationofR2¼0:5,however,removingthree entificobjectives.Inthesecaseswebelieveharvestedmaterialwas
eithermixedfrommultipletreesinthecaseofthelargeoutlier,or
significant outliers resulted in R2¼0:75 with 28 samples, which
erroneously taken from adjacent trees in the case of the smaller
is similar to first year. The data, outliers and linear relationship
outliers.
are shown in Fig. 11(a). It was considered reasonable to reject
ThelidarvolumesandyieldsforNonpareiltreesforbothyears
thethreepointsduetothecombinationofthefollowingreasons:
areshowntogetherinFig.11(b),withatotalof11and31samples
Theyare allat least twiceas far from the lineof best fit than all
respectively.Datafrombothyearsareonthesametrendline.Ifwe
other points. They represent the two highest recorded absolute
rejectthesamethreeoutliersasabove,asinglelineofbestfitmod-
harvestweights(frompost-harvestweighing)andthesinglelow-
elsbothyearswell,withR2¼0:77for39samples,whichissignif-
estmeasuredcanopyvolumerespectively(fromoursystem).The
icant.Thisequatestoastandarddeviationof6.26kgpertree.Over
latterwasconsistentlymeasuredtobeasmall,replantedtreeboth
asimilarorchardblockwith580trees,theexpectedstandarddevi-
years,whichisunlikelytoyield50kginthesecondyear.Consid-
ationoftheestimatedsum-totalyieldisapproximately150kg,but
eringonlytheharvestdata(and notthemeasurementsfromour
significantly larger experiments over multiple orchard blocks
system),for these outliers to be true would implythat two trees
would be required to characterise the whole-block error accu-
rately. Furthermore, whole-block sample weights might provide
amethodtocalibrateandconstraintheerror,butthisremainsas
9 Increasingthesamplesizewaschallengingduetothelabourinvolvedinselective
harvestingatthecommercialorchard. futureworktovalidate.Thefigurealsoconnectsthe11Nonpareil
92 J.P.Underwoodetal./ComputersandElectronicsinAgriculture130(2016)83–96
Harvest Weight vs Lidar Canopy Volume (R 2=0.75) Canopy Volume vs Harvest Yield Over Two Years (R2=0.77)
70 70
60 60
3m) 3m)
me ( 50 me ( 50
Canopy Volu 3400 Canopy Volu 3400 YYeeaarr 12:: 22001145
Lidar 20 Lidar 20
10 10
0 0
0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90
Harvest Weight (kg) Harvest Weight (kg)
Fig.11. Lidarcanopyvolumeestimatesandyieldfor31Nonpareiltrees(a)in2014and(b)2014and2015together,with11treesmeasuredbothyearsconnectedwithgrey
lines.Threesamplesareassumedtobeoutliers,circledinred,whicharenotincludedinthelinefitorR2calculation.(Forinterpretationofthereferencestocolourinthis
figurelegend,thereaderisreferredtothewebversionofthisarticle.)
Fig.12. Amapshowingthepredictedharvestyieldinthesecondyear,fromthelidarcanopyvolumesmeasuredatflowering.Themeanpredictedyieldis46.4kgpertree.
Blackthroughred,yellow,towhiterepresentstherangefrom30to60kg.Theredareatotheleftofthefigureshowsaregionofconsistentlyloweryield,whichcorresponds
tothegrowersassessmentthatlessyieldisproducednearthelargeaccessroad.Thesedataallowtheeffecttobequantified.(Forinterpretationofthereferencestocolourin
thisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
trees measured both years, showing that the volume and yield the southern (left) region near the access road. This was already
both tended to increase from 2014 to 2015. The data show that knowntothefarmmanager,however,oursystemisabletoquan-
lidar-canopy volume can explain most of the variation in yield, tifyanddisplaythemagnitudeoftheeffect,whichmaybeauseful
bothwithinasingleseasonwhenconsideringdifferenttreesand considerationformanagementdecisions.
alsofromoneyeartothefollowingyear.Thetrendsuggeststhat
thelidarvolumemaybeagoodpredictororproxyforyield,even
3.3.Flowerdensityconsistency
whenmeasuredearlyintheseasonduringflowering,butalarger
samplesizecoveringadditionalyearsanddifferentorchardswould
The consistency of the image based flower density estimates
beprudenttodrawbroaderconclusions.Inparticular,otherfactors
was assessed by comparing the measurements of the same trees
such as variable pollination, weather conditions or disease that
on subsequent days at peak bloom, in both years. The peak time
werenotcontrolledforinourexperimentscouldalterthestrength
wasestimatedvisuallybythefarmmanager.Inthefirstyear,rel-
oftherelationships.
The linear model from both years was applied to the lidar ativelypooragreementandrepeatabilitywasobserved(R2¼0:61,
canopyvolumestakenduringfloweringinthesecondyear,topro- CV=17.9%),duetodifficultieswiththenaturalillumination.Five
duceamapofthepredictedharvestyieldforallscannedNonpareil of the ten rows were scanned when the angle to the sun caused
trees,showninFig.12.Themodeliscalibratedbythelinearfitto lens flaring and partial saturation, which degraded the quality of
selectively harvested individual trees from both years, and then theimages.Thefiveunaffectedrows(roughly290trees)hadbetter
usedtoestimatetheyieldinkilogramsforeverytreeinthesecond consistency (R2¼0:79, CV=11.0%). The illumination conditions
year. The relatively small number of individual sample trees are wereimprovedinthesecondyear,byscanningtheeastandwest
usedtofitthemodel,buttheneverytreecanbemoreaccurately faces, in the morning and afternoon with the sun always behind
estimated,comparedtoaprocessofinterpolationorextrapolation thecamera.Alltenrows(roughly580trees)wereusedinthesec-
fromthesparsesamples. ond year, and consistency improved (R2¼0:85, CV=9.1%). The
Qualitatively, the map shows good alignment between flower density repeatability is lower than the lidar volume
replanted trees and the gaps in the canopy seen in the satellite (CV=2.4%),whichisduetothedynamicsoftheflowers:onadaily
photo.Furthermore,thereisavisibleregionofdepressedyieldin time-scale, flowers can increasedue to continuous bloomingand