void edgeDetection(){
	Mat src;
	char fname[MAX_PATH];

	while (openFileDlg(fname)){
		src = imread(fname, CV_LOAD_IMAGE_GRAYSCALE);

		Mat gauss = gaussian2D(src, 5, 0.8);

		Mat gradient = Mat(gauss.rows, gauss.cols, CV_32FC1);
		Mat fi = Mat(gauss.rows, gauss.cols, CV_32FC1);
		Mat gradientShow = Mat(gauss.rows, gauss.cols, CV_8UC1);
		Mat nonMaxima = Mat(gauss.rows, gauss.cols, CV_32FC1);

		float gx, gy, G, FI;
		int area;

		for (int i = 1; i < gauss.rows - 1; i++)
			for (int j = 1; j < gauss.cols - 1; j++){
				gx = 0;
				gy = 0;

				gx += -gauss.at<uchar>(i - 1, j - 1);
				gx += 0;
				gx += gauss.at<uchar>(i - 1, j + 1);
				gx += - 2 * gauss.at<uchar>(i, j - 1);
				gx += 0;
				gx += 2 * gauss.at<uchar>(i, j + 1);
				gx += -gauss.at<uchar>(i + 1, j - 1);
				gx += 0;
				gx += gauss.at<uchar>(i + 1, j + 1);

				gy += gauss.at<uchar>(i - 1, j - 1);
				gy += 2 * gauss.at<uchar>(i - 1, j);
				gy += gauss.at<uchar>(i - 1, j + 1);
				gy += 0;
				gy += 0;
				gy += 0;
				gy += -gauss.at<uchar>(i + 1, j - 1);
				gy += -2*gauss.at<uchar>(i + 1, j);
				gy += -gauss.at<uchar>(i + 1, j + 1);

				G = sqrt(gx*gx + gy*gy);
				G /= (4 * sqrt(2));
				FI = atan2(gy, gx);


				gradient.at<float>(i, j) = G;
				fi.at<float>(i, j) = FI;
				gradientShow.at<uchar>(i, j) = (uchar)G;
			}
		nonMaxima = gradient.clone();
		for (int i = 1; i < fi.rows - 1; i++)
			for (int j = 1; j < fi.cols - 1; j++){
				area = circleZone(fi.at<float>(i, j));
				switch(area){
				case 0:
					if (gradient.at<float>(i, j - 1) > gradient.at<float>(i, j) || gradient.at<float>(i, j + 1) > gradient.at<float>(i, j))
						nonMaxima.at<float>(i, j) = 0;
					break;
				case 1:
					if (gradient.at<float>(i - 1, j + 1) > gradient.at<float>(i, j) || gradient.at<float>(i + 1, j - 1) > gradient.at<float>(i, j))
						nonMaxima.at<float>(i, j) = 0;
					break;
				case 2:
					if (gradient.at<float>(i - 1, j) > gradient.at<float>(i, j) || gradient.at<float>(i + 1, j) > gradient.at<float>(i, j))
						nonMaxima.at<float>(i, j) = 0;
					break;
				case 3:
					if (gradient.at<float>(i - 1, j - 1) > gradient.at<float>(i, j) || gradient.at<float>(i + 1, j + 1) > gradient.at<float>(i, j))
						nonMaxima.at<float>(i, j) = 0;
					break;
				default:
					break;
				}


			}

		Mat nonMaximaShow = Mat(nonMaxima.rows, nonMaxima.cols, CV_8UC1, Scalar(0));
		for (int i = 1; i < fi.rows - 1; i++)
			for (int j = 1; j < fi.cols - 1; j++){
				nonMaximaShow.at<uchar>(i, j) = (uchar)nonMaxima.at<float>(i, j);
			}

		int noEdgePixels = 0;

		for (int i = 1; i < nonMaximaShow.rows - 1; i++)
			for (int j = 1; j < nonMaximaShow.cols - 1; j++){
				if (nonMaximaShow.at<uchar>(i, j) > 0)
					noEdgePixels++;
			}
		noEdgePixels *= 0.1;
		int* histogram = generate_histogram(nonMaximaShow);
		int  noEdgePixelsSum = 0;

		int tHigh = 255;
		for (tHigh = 255; noEdgePixelsSum < noEdgePixels; noEdgePixelsSum += histogram[tHigh], tHigh--);
		int tLow = 0.4 * tHigh;

		for (int i = 1; i < nonMaximaShow.rows - 1; i++)
			for (int j = 1; j < nonMaximaShow.cols - 1; j++){
				if (nonMaximaShow.at<uchar>(i, j) >= tHigh)
					nonMaximaShow.at<uchar>(i, j) = 255;
				else if (nonMaximaShow.at<uchar>(i, j) <= tHigh && nonMaximaShow.at<uchar>(i, j) >= tLow)
					nonMaximaShow.at<uchar>(i, j) = 127;
				else
					nonMaximaShow.at<uchar>(i, j) = 0;
			}

		//Edge linking
		int offi[9] = { -1, -1, -1, 0, 0, 0, 1, 1, 1 };
		int offj[9] = { -1, 0, 1, -1, 0, 1, -1, 0, 1 };

		for (int i = 1; i < nonMaximaShow.rows - 1; i++){
			for (int j = 1; j < nonMaximaShow.cols - 1; j++){
				if (nonMaximaShow.at<uchar>(i, j) == 255){
					queue<Point2i> Q;
					Q.push(Point2i(j, i));

					while (!Q.empty()){
						Point2i q = Q.front();
						Q.pop();
						for (int k = 0; k < 9; k++){
							if (nonMaximaShow.at<uchar>(q.y + offi[k], q.x + offj[k]) == 255){
								labels.at<int>(q.y + offi[k], q.x + offj[k]) = label;
								Q.push(Point2i(q.x + offj[k], q.y + offi[k]));
							}
						}
					}
				}
			}
		}

		imshow("Step2", nonMaximaShow);
		waitKey();

	}
}