本发明涉及林木植株的调查与研究技术领域，具体涉及一种无人机图像的林木植株并行提取方法。
背景技术：
无人机被动光学遥感在时空分辨率、可获取性、成本方面较卫星被动光学遥感具有较大优势，当前越来越多的被应用于林木植株的调查与研究中。对于林木(主要是果园、人工林等)植株提取而言，主要是利用数字图像处理技术，从无人机拍摄的图像中识别出目标植株个体并进行统计或分析。
然而，目前的数字图像处理技术的精确度普遍较低，不同冠层的林木难以同时识别。同时，由于当前无人机图像分辨率普遍较高，使得数字图像处理需要进行大量的运算，导致算法耗时较长。
技术实现要素：
针对现有技术的不足，本发明的目的在于提供一种无人机图像的林木植株并行提取方法，以提高图像提取的精度和速度。
为了实现上述目的，本发明采取的技术方案是：
一种无人机图像的林木植株并行提取方法，包括步骤：
获取无人机拍摄的林区图像；
采用GPU并行处理的方式，对林区图像进行基于尺度空间技术的二进制大型对象检测，得到初始检测点；
对初始检测点进行筛选，将不含绿色的检测点删除；
将筛选后的检测点作为提取结果。
与现有技术相比，本发明的有益效果在于：
本发明采用尺度空间技术对无人机获取的林区图像进行检测，可以在各种尺度下识别不同冠层大小的林木，由于加入了对象分析步骤，精度较高。同时，对关键步骤采用GPU并行处理的方式，使得本提取方法具有较高的效率和速度，为林木植株的调查与研究提供了支持。
附图说明
图1为本发明无人机图像的林木植株并行提取方法的流程示意图。
具体实施方式
下面结合具体实施方式对本发明作进一步的说明。
本发明无人机图像的林木植株并行提取方法，如图1所示，包括步骤：
步骤s101、获取无人机拍摄的林区图像；
步骤s102、采用GPU并行处理的方式，对林区图像进行基于尺度空间技术的二进制大型对象检测，得到检测点；
步骤s103、对初始检测点进行筛选，将不含绿色的检测点删除；
步骤s104、将筛选后的检测点作为提取结果。
利用PyCUDA平台，实现无人机图像的林木植株并行提取算法。CUDA(Compute Unified Device Architecture)技术为NVIDIA公司开发的类C语言GPU编程平台。而PyCUDA是CUDA的Python语言封装，该封装提供对象自动清理功能，并通过结合Numpy、Scipy和GDAL等Python科学计算库，能够方便的开发复杂应用。
上述算法可以分为：1)输入；2)基于尺度空间技术的二进制大型对象(Binary Large Objects，BLOBS)检测；3)对象分析；4)输出，共四个部分。其中第3部分为图像处理，能够进行GPU并行处理，并且可以取得效率提升，是本发明的核心部分，采用并行算法编程，其余部分均采用串行编程。
其中，BLOBS检测部分包含3个核心步骤：
1)原始图像的多尺度LoG(Laplacian of Gaussian，LoG，拉普拉斯-高斯)滤波；
2)LoG滤波尺度空间图像的局部最大值滤波；
3)检测点的基于面积的组合筛选。
下面详细介绍上述步骤。
算法总共包含5个参数，分别为m_grey,min_sigma,max_sigma,num_sigma,threshold,overlap。(注：Sigma代表高斯滤波的尺度，其值越高模糊程度越高)其中，m_grey是输入的灰度图像，min_sigma,max_sigma,num_sigma分别对应拉普拉斯-高斯滤波的最小尺度、最大尺度和尺度数量，threshold代表BLOB识别灰度阈值，overlap代表BLOB筛选面积重叠阈值。
技术内容将按照算法处理步骤分别阐述：
(1)输入：通过GDAL读入无人机可见光RGB(Red,Green,Blue，红绿蓝)图像。利用公式1将RGB波段转换为灰度图像，Y表示灰度值。
Y＝0.2125*R+0.7154*G+0.0721*B (公式1)
(2)基于尺度空间技术的二进制大型对象(Binary Large Objects，BLOBS)检测：
准备步骤：根据参数min_sigma,max_sigma,num_sigma，分别计算num_sigma个等距的sigma值，记为sigma_list
1)原始图像的多尺度LoG滤波：针对原始图像m_grey，分别采用sigma_list中的各个sigma进行LoG滤波，得到LoG滤波尺度空间图像。其中，LoG滤波的公式为：
其中，x,y分别为图像中的x,y坐标，σ(sigma)为高斯滤波尺度。
实现流程为：
由于处理无人机图像需要较大的sigma(≤50)，而当前CUDA设备允许的一个块(Block)中的线程(Thread)数不能超过1024个(32*32)，当sigma为8.0是窗口大小已达到32。因此，二维滤波器无法满足需求，必须采用两个一维LoG滤波器进行组合，以便sigma较大时进行LoG滤波。组合方法为，针对某一sigma：
①对原始图像进行列方向LoG滤波，然后进行行方向高斯滤波，得到Mlogx；
②对原始图像进行行方向LoG滤波，然后进行列方向高斯滤波，得到Mlogy；
③Mlog＝(Mlogx+Mlogy)*(-sigma2)。
最终针对sigma_list中的每个sigma进行滤波，并合并为三维图像数组。
注：Kernel函数采用常数内存(Constant Memory)存储滤波器参数；而对于每行或列采用一维数组索引共享内存(Shared Memory)，用以加快运算速度，即每个Block内的所有Thread(1024个)共用一份缓存的图像数据。由于滤波器尺寸较大，为避免边缘数据浪费，边界处理方式为反射(Reflect)，即M(x-i)＝M(x+i-1)。第③步由于不涉及滤波，直接在主显示内存(Global Memory)操作。
2)LoG滤波尺度空间图像的局部最大值滤波
实现流程为：
①窗口尺寸为3(半径为1)个像素，分别求取三维LoG滤波图像组在x,y,z方向的最大值，并生成最大值立方体；
②对比大值立方体和LoG立方体，对于满足公式的像素作为潜在的BLOBS。
ImageLoG(x,y,z)＝Imagemax(x,y,z)并且ImageLoG(x,y,z)>Threshold (公式2)
③输出所有潜在BLOBS的x,y坐标，以及z轴对应的sigma值
注：对于每行或列采用一维数组索引共享内存(Shared Memory)，用以加快运算速度，即每个Block内的所有Thread(1024个)共用一份缓存的图像数据。而对于z轴，由于数量小于1024，与x轴共同分配1024个Threads，采用二维数组索引共享内存。为避免滤波器在边缘时窗口像元数量减少导致过多的像素被选取，边界处理方式为常数(Constant)，即M(x-i)＝0。
3)检测点的基于面积的组合筛选
输入潜在的BLOBS，为Python语言的List格式，每个单元包含x,y,sigma三元素。本步将对所有BLOBS进行两两组合，对比BLOBS的面积，如果某两个BLOB的重叠度小于阈值overlap，则删除面积较小的BLOB。本步骤利用GPU实现的难点在于，由于BLOBS长度未知，难以采用二维数组用Thread索引组合中的两个BLOB元素。本发明采用一维数组索引组合序号，然后通过组合序号反推BLOB元素。原理为代码1：
实现流程为：
①根据BLOBS数量NBLOBS，利用公式3在CPU端计算两两组合的总数Ncombination
Ncombination＝NBLOBS*(NBLOBS-1)/2 (公式3)
②启动一维数组，为每个Block,Thread计算其对应的组合一维序号，采用迭代法计算其对应的两个BLOBS序号；
③利用公式4计算重叠度OL：
R1＝sigma1*21/2,R2＝sigma2*21/2
Distance＝((x1–x2)2+(y1–y2)2)1/2
如果Distance>R1+R2
OL＝0
如果Distance≤|R1-R2|
OL＝1
否则：
Ratio1＝(Distance2+R12–R22)/(2.0*Distance*R1)
且Ratio1＝max(Ratio1,-1.0)，Ratio1＝min(Ratio1,1.0)
Acos1＝arccos(Ratio1)
Ratio2＝(Distance2+R22–R12)/(2.0*Distance*R2)
且Ratio2＝max(Ratio2,-1.0)，Ratio2＝min(Ratio2,1.0)
Acos2＝arccos(Ratio2)
A＝-Distance+R2+R1
B＝Distance-R2+R1
C＝Distance+R2-R1
D＝Distance+R2+R1
Area＝R12*Acos1+R22*Acos2-0.5*(|A*B*C*D|)1/2
Amin＝min(R1,R2)
OL＝Area/(π*Amin2) (公式4)
式中，sigma1和sigma2分别表示高斯滤波尺度1和高斯滤波尺度2，x1、y1表示高斯滤波尺度1下，潜在检测点1的坐标，x2、y2表示高斯滤波尺度2下，潜在检测点2的坐标，Distance表示欧氏距离；R1、R2、Ratio1、Ratio2、Acos1、Acos2、A、B、C、D、Area、Amin均表示中间变量。
④在CPU端，删除所有OL值小于overlap阈值的BLOB，剩余的BLOBS即为提取的BLOBS，接下来通过后续步骤分析其颜色，进一步筛选。
注：由于GPU芯片的对称性，只有当分配的Blocks*Threads值为2的整数次方时会发挥最佳性能，而GTX960m显卡允许的Blocks*Threads最大值为2097152(221)，因此本发明将该值作为Kernel启动参数。
(3)对象分析：将原始图像从RGB色彩空间转换为CIE-Lab色彩空间，转换方式见公式9。根据上一步提取的BLOBS位置和面积，计算所有BLOB的a通道均值，删除值高于0的BLOBS(a值大于0表示不含绿色)，剩余点为检测结果
X＝R*0.4124+G*0.3576+B*0.1805 (公式9)
Y＝R*0.2126+G*0.7152+B*0.0722
Z＝R*0.0193+G*0.1192+B*0.9505
L＝116*f(Y)–16
a＝500*(f(X)-f(Y))
b＝200*(f(Y)-f(Z))
其中，当x>0.008856时，f(x)＝x^(1/3)，x表示X、Y、Z中的任意一个，当x<＝0.008856时，f(x)＝(7.787*x)+(16/116)。
(4)输出：输出为BLOBS列表，其中包含每个BLOBS的x,y坐标及其半径(注：半径为sigma)。
上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。