/*
函数功能：使用阈值分割，双线性插值，分割亚像素边缘
函数参数：输入的灰度图，灰度阈值
函数返回值：无
*/
void threshold_sub_pixel(Mat src, int th)
{
	int border_r = 0,border_c=0;  //像素级别边缘点个数，先留着
	int flag = 0;
	int i, j;
	double f = 0.0, x=0.0, y=0.0; 
	double sample_num = 0.05; // 采样精度，一般为0.05或者0.04
	double errors = 0.2; //允许误差，这个值的选择有待考究
	double percent = 0.9;  //允许有（1-percent）*100%的像素边缘丢失
	int Q11, Q22, Q12, Q21, max, min;

	//FILE* fp = NULL;
	//fp = fopen("C:\\Users\\13553474338\\Desktop\\graduate_design\\sub_pixel.xls", "w");
	while ((double)border_c / (border_r + 1) < percent) //需要满足一定丢失率
	{
		for (i = 0; i < src.rows - 1; i++)
		{
			for (j = 0; j < src.cols - 1; j++)
			{
				Q12 = src.at <uchar>(i, j);
				Q22 = src.at <uchar>(i, j + 1);
				Q11 = src.at <uchar>(i + 1, j);
				Q21 = src.at <uchar>(i + 1, j + 1);
				max = my_max(Q11, Q12, Q21, Q22);
				min = my_min(Q11, Q12, Q21, Q22);
				//如果阈值落在某个2*2的区域内，那么阈值必须小于四个灰度值中最大的，并且大于四个灰度值最小的
				if (th< max + 1 && th > min - 1)
				{
					border_r++;
					flag = 0;
					for (x = 0; x < 1 + sample_num && flag == 0; x += sample_num)
					{
						for (y = 0; y < 1 + sample_num && flag == 0; y += sample_num)
						{
							f = Q11 * (1 - x) * (1-y) + Q21 * x * (1-y) + Q12 * (1 - x) * y  + Q22 * x * y;
							if (abs(f - th) < errors)
							{
								//std::cout << "像素边缘坐标为：(" << i << "，" << j << ")，" << "亚像素边缘坐标为：(" << i + x << "，" << j + y << ")" << endl;
								//fprintf(fp, "%f\t%f\t%d\n", j+x+1, i+1-y+1, src.at <uchar>(i, j)); //保存数据
								border_c++;
								flag = 1;
							}
							if (flag == 1)
							{
								continue;
							}
						}
						if (flag == 1)
						{
							continue;
						}
					}
				}
			}
		}
		std::cout << endl << "像素边缘点数：" << border_r << "，亚像素边缘点数为：" << border_c << ", 可以找到亚像素边缘点的概率为" << (double)border_c * 100 / border_r << "%" << endl;
		errors += errors;
	}
	//fclose(fp);
	return;
}