首先来看原图像
原图像中光干扰较大，因此我们后面的处理会使用下图红框这样只有杆影的局部区域进行计算。
图像切割拆分
视频的图像帧数为25帧长度为从8:54:06到9:34:46共40min,把图像每隔3min切割一张(3min=25*180=4500帧)并保存
clcclearobj = VideoReader('Appendix4.avi')numFrames = obj.NumberOfFramesfor k = 1 :4500: numFrames frame = read(obj,k); imshow(frame); %imwrite(frame,strcat(num2str(k),'.jpg'),'jpg');end 图像处理
Canny边缘检测
Canny边缘检测算法的处理流程
Canny边缘检测算法可以分为以下5个步骤：
- 使用高斯滤波器，以平滑图像，滤除噪声。
- 计算图像中每个像素点的梯度强度和方向。
- 应用非极大值（Non-Maximum Suppression）抑制，以消除边缘检测带来的杂散响应。
- 应用双阈值（Double-Threshold）检测来确定真实的和潜在的边缘。
- 通过抑制孤立的弱边缘最终完成边缘检测。 下面详细介绍每一步的实现思路。
高斯平滑滤波
为了尽可能减少噪声对边缘检测结果的影响，所以必须滤除噪声以防止由噪声引起的错误检测。为了平滑图像，使用高斯滤波器与图像进行卷积，该步骤将平滑图像，以减少边缘检测器上明显的噪声影响。大小为(2k+1)x(2k+1)的高斯滤波器核的生成方程式由下式给出： 
下面是一个sigma = 1.4，尺寸为3x3的高斯卷积核的例子（需要注意归一化）：
若图像中一个3x3的窗口为A，要滤波的像素点为e，则经过高斯滤波之后，像素点e的亮度值为：
其中*为卷积符号，sum表示矩阵中所有元素相加求和。 高斯卷积核大小的选择将影响Canny检测器的性能。尺寸越大，检测器对噪声的敏感度越低，但是边缘检测的定位误差也将略有增加。 计算梯度强度和方向
图像中的边缘可以指向各个方向，因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。边缘检测的算子（如Roberts，Prewitt，Sobel等）返回水平Gx和垂直Gy方向的一阶导数值，由此便可以确定像素点的梯度G和方向theta。
其中G为梯度强度， theta表示梯度方向，arctan为反正切函数。下面以Sobel算子为例讲述如何计算梯度强度和方向。
x和y方向的Sobel算子分别为：
其中Sx表示x方向的Sobel算子，用于检测y方向的边缘； Sy表示y方向的Sobel算子，用于检测x方向的边缘（边缘方向和梯度方向垂直）。在直角坐标系中，Sobel算子的方向如下图所示。
若图像中一个3x3的窗口为A，要计算梯度的像素点为e，则和Sobel算子进行卷积之后，像素点e在x和y方向的梯度值分别为：
其中*为卷积符号，sum表示矩阵中所有元素相加求和。根据公式（3-2）便可以计算出像素点e的梯度和方向。
非极大值抑制
非极大值抑制是一种边缘稀疏技术，非极大值抑制的作用在于“瘦”边。对图像进行梯度计算后，仅仅基于梯度值提取的边缘仍然很模糊。对于标准3，对边缘有且应当只有一个准确的响应。而非极大值抑制则可以帮助将局部最大值之外的所有梯度值抑制为0，对梯度图像中每个像素进行非极大值抑制的算法是：
- 将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。
- 如果当前像素的梯度强度与另外两个像素相比最大，则该像素点保留为边缘点，否则该像素点将被抑制。 通常为了更加精确的计算，在跨越梯度方向的两个相邻像素之间使用线性插值来得到要比较的像素梯度，现举例如下：
如上图所示,将梯度分为8个方向，分别为E、NE、N、NW、W、SW、S、SE，其中0代表00~45o,1代表450~90o，2代表-900~-45o，3代表-450~0o。像素点P的梯度方向为theta，则像素点P1和P2的梯度线性插值为：
非极大值抑制的伪代码描写如下：
梯度方向的计算要和梯度算子的选取保持一致。
遍历当前目录所有截图并且读取图像并且提取兴趣区域(裁剪出只有杆子影子的图像区域)，对图像进行灰度化处理，方便转化为二值化图像，转化为二值图像以后使用Canny边缘检测算法，提取图像边缘，过滤噪声图像区域只保留影子变化区域。
经过上述步骤后从右至左遍历图像，寻找最右的像素坐标点即为影子顶点。过滤噪声前图像:
过滤噪声后图像:
clcclearcount=1for i=1:4500:58501 %disp(strcat(num2str(i),'.jpg')); BW=imread(strcat(num2str(i),'.jpg')); BW=imcrop(BW,[922 850 800 90]); BW=rgb2gray(BW); thresh=[0.1,0.17]; sigma=1.4;%定义高斯参数 %Canny边缘检测 f = edge(double(BW),'canny',thresh,sigma); for i=50:1:91 for j=520:1:900 f(i,j)=0; end end %figure(1),imshow(imcrop(f,[520 40 1000 80]),[]); [f,x,y]=find_pwite(f); dis_x(1,count)=x; dis_y(1,count)=y; count=count+1; %imwrite(f,strcat(num2str(count),'.jpeg'),'jpg'); figure(i),imshow(f,[]); %[x,y,button] = ginputend
绘制时间变化图像
for i=1:length(dis_x); h(1,i)=power(power(dis_x(1,i),2)+power(dis_y(1,i),2),1/2);endt=8.9:0.05:9.55;plot(t,h,'r');
