极地遥感学习笔记（一）——冰间水道

本文学习内容参考文献：

[1]王玉,马玉贤,陈元,陈伟斌.基于卫星图像的北极冰间水道形态学特征分析[J].海洋技术学报,2019,38(02):8-13.
[2] Cunningham G F, Kwok R, Banfield J. Ice lead orientation characteristics in the winter Beaufort Sea [C]// International Geoscience & Remote Sensing Symposium, IEEE, 1994
[3]Tschudi M A, Curry J A, Maslanik J A. Characterization of springtime leads in the Beaufort/Chukchi Seas from airborne and satellite
observations during FIRE/SHEBA[J]. Journal of Geophysical Research Oceans, 2002, 107(C10):SHE 9-1-SHE 9-14
[4]Banfield J. Skeletal modeling of ice leads[J]. IEEE Transactions on Geoscience & Remote Sensing, 1992, 30(5):918-923

一、提取冰间水道的骨架

1.1 前提：

①首先对研究区的遥感影像进行几何校正、云掩膜等预处理（由于极地地区获取的光学影像存在着大量的云遮盖，要想将雾与云区分出来必须进行云掩膜）
②运用劈窗法计算冰面温度。统计直方图确定阈值，通过阈值分割将影像处理成一副二值图象，经过 Hilditch算法提取出冰间水道的骨架。

1.2 方法 ：

Hilditch算法

1.2.1 原理：

算法主要步骤为遍历图像像素。在一次遍历中，对每个像素，若同时满足6个条件，则将其标记为待移除，一次遍历后若存在待移除像素，则将本次遍历所有待移除像素值更改。然后重新开始遍历直到不存在待移除像素，遍历结束。

1.2.2 六个判定条件：

①当前像素点P值为1，即不是背景。
②P的上、下、左、右邻居至少有一个值为0。
③P的八邻居中至少有2个为1。
④P的八邻居中至少有一个邻居值为1而未被标记为待移除。
⑤计算P的连接数（两种方案选一种即可）:

• 4连通连接数计算公式：

• 8连通连接数计算公式：
  
  ⑥若x3和x5已标记为待移除，分别用x3=0和x5=0重新计算连接数且要等于1。
  Hilditch参考原文链接: link.

1.3 opencv代码

#include <iostream>
#include <opencv2\opencv.hpp>
 
using namespace std;
using namespace cv;
int Judge_Condition1(Mat &src, int i, int j);
int Judge_Condition2(Mat &src, int i, int j);
int Judge_Condition3(Mat &src, int i, int j);
int Judge_Condition4(Mat &mask, Mat &src, int i, int j);
int Judge_Condition5(Mat &mask, Mat &src, int i, int j);
struct Direct
{
	int x;
	int y;
};  //方向结构体
	//定义n0...n7八个方向
Direct N[8] = { { -1,0 },{ -1,1 },{ 0,1 },{ 1,1 },{ 1,0 },{ 1,-1 },{ 0,-1 },{ -1,-1 } };
int main(void)
{
	//读取图片
	Mat src = imread("K:/line.jpg");
	Mat image;
	//转为灰度图像
	cvtColor(src, image, CV_BGR2GRAY);
	int rows = image.rows;
	int cols = image.cols;
	
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			if ((int)image.at<uchar>(i, j) > 200)
			{
				image.at<uchar>(i, j) = 255;
			}
			else
			{
				image.at<uchar>(i, j) = 0;
			}
		}
	}
	//复制图像
	Mat dst_temp = image;
	Mat mask=image.clone();
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			mask.at<uchar>(i, j) = 0;
		}
	}
	namedWindow("src", 1);
	imshow("src", image);
	bool flag = true;
	while(flag)
	{
		flag = false;
		int count = 0;
		for (int i = 1; i < rows - 1; i++)
		{
			for (int j = 1; j < cols - 1; j++)
			{
				if (image.at<uchar>(i, j) == 0)
					continue;
				if (Judge_Condition1(image, i, j) == 1)
					continue;
				if (Judge_Condition2(image, i, j) == 0)
					continue;
				if (Judge_Condition3(image, i, j) != 1)
					continue;
				if (Judge_Condition4(mask, image, i, j) == 1)
					continue;
				if (Judge_Condition5(mask, image, i, j) == 1)
					continue;
				mask.at<uchar>(i, j) = 1;
				count++;
			}
		}
		if (count)
		{
			for (int i = 0; i < mask.rows; i++)
			{
				for (int j = 0; j < mask.cols; j++)
				{
					if (mask.at<uchar>(i, j) == 1)
						image.at<uchar>(i, j) = 0;
				}
			}
			flag = true;
			/*************一次遍历后将标记复位*******************/
			for (int i = 0; i < image.rows; i++)
			{
				for (int j = 0; j < image.cols; j++)
				{
					mask.at<uchar>(i, j) = 0;
				}
			}
		}
	}
	
	
	namedWindow("dst", 1);
	imshow("dst", image);
	waitKey(0);
	return 0;
}
 
int Judge_Condition1(Mat &image, int i, int j)//四邻域点不全为1
{
	if (((int)image.at<uchar>(i + N[0].x, j + N[0].y) == 0 )&&
		((int)image.at<uchar>(i + N[2].x, j + N[2].y) == 0 )&&
		((int)image.at<uchar>(i + N[4].x, j + N[4].y) == 0 )&&
		((int)image.at<uchar>(i + N[6].x, j + N[6].y) == 0))
	{
		return 1;
	}
	else
	{
		return 0;
	}
}
 
int Judge_Condition2(Mat &image, int i, int j)
{
	int sum = 0;
	sum = (int)image.at<uchar>(i + N[0].x, j + N[0].y) +
		(int)image.at<uchar>(i + N[1].x, j + N[1].y) +
		(int)image.at<uchar>(i + N[2].x, j + N[2].y) +
		(int)image.at<uchar>(i + N[3].x, j + N[3].y) +
		(int)image.at<uchar>(i + N[4].x, j + N[4].y) +
		(int)image.at<uchar>(i + N[5].x, j + N[5].y) +
		(int)image.at<uchar>(i + N[6].x, j + N[6].y) +
		(int)image.at<uchar>(i + N[7].x, j + N[7].y);
	if (sum >= 510)
	{
		return 1;
	}
	else
	{
		return 0;
	}
}
 
int Judge_Condition3(Mat &image, int i, int j)
{
	int neighbor[8] = { 0 };
	for (int k = 0; k < 8; k++)
	{
		neighbor[k] = (int)image.at<uchar>(i + N[k].x, j + N[k].y) == 0 ? 1 : 0;
	}
	int count = neighbor[0] - (neighbor[0] & neighbor[1] & neighbor[2]);
	count += neighbor[2] - (neighbor[2] & neighbor[3] & neighbor[4]);
	count += neighbor[4] - (neighbor[4] & neighbor[5] & neighbor[6]);
	count += neighbor[6] - (neighbor[6] & neighbor[7] & neighbor[0]);
	return count;
 
}
 
int Judge_Condition4(Mat &mask, Mat &image, int i, int j)
{
	if (mask.at<uchar>(i - 1, j) == 1)
	{
		image.at<uchar>(i - 1, j) = 0;
		if (Judge_Condition3(image, i, j) != 1)
		{
			image.at<uchar>(i - 1, j) = 255;
			return 1;
		}
	}
}
 
int Judge_Condition5(Mat &mask, Mat &image, int i, int j)
{
	if (mask.at<uchar>(i, j - 1) == 1)
	{
		image.at<uchar>(i, j - 1) = 0;
		if (Judge_Condition3(image, i, j) != 1)
		{
			image.at<uchar>(i, j - 1) = 255;
			return 1;
		}
	}
}

二、计算冰间水道的倾角

2.1 方法：线段拟合、霍夫变换

2.1.1 原理：

继得到的冰间水道中心骨架利用线段进行拟合，即将整个骨架分割成若干近似于直线的线段。利用霍夫变换计算线段的斜率即冰间水道的倾角。
参考原文(opencv计算角度)链接:link

三、计算冰间水道的距离

方法：

同样利用线段对冰间水道的骨架进行拟合，确定对应线段的中心点即为冰
间水道的中心点位置。计算每个冰间水道中心点与其他冰间水道中心点的距离，寻找最短距离，即为该冰间水道与其相邻的冰间水道之间的距离。计算所有相邻冰间水道之间的距离的算数平均值作为图像中相邻冰间水道的平均距离。

四、计算冰间水道的长度

方法：

在计算冰间水道的长度时，Cunningham [2]提供了 3 种冰间水道长度的测量
主对角线长度、冰间水道骨架长度以及冰间水道总长度。文中采用主对角线长度为连接骨架最远的两个端点之间的距离，冰间水道骨架长度为沿着骨架连接骨架端点的距离，冰间水道总长度为骨架长度加上分支的长度。由于冰间水道的分支对冰间水道的影响较小，因此采用主对角线长度来表示冰间水道的长度。

参考代码链接：link