이미지처리 잡다한것들

#include "utils.h"

Mat Histogram_equalization(Mat src)
{
	/*
	* 1. 입력 영상의 히스토그램 획득
	* 2. 히스토그램을 정규화(확률화)
	* 3. 변환 테이블을 만들기 위해, 정규화된 히스토그램으로부터 각 밝기별로 CDF를 계산한다.
    * 4. STEP 3을 통해 얻은 변환 테이블을 기반으로 r -> s로의 갱신을 수행한다.
	*/
	Mat dst = src; 
	int min = 2100000000, max = 0;
	int totalpixel = src.rows * src.cols; 

	//히스토그램 구함
	int hist[256] = { 0, };
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			hist[dst.at<uchar>(i, j)] ++;
			if (src.at<uchar>(i, j) < min) min = src.at<uchar>(i, j);
			if (src.at<uchar>(i, j) > max) max = src.at<uchar>(i, j);
		}
	}

	//평활화, 빈도값 구하기
	double equalhist[256] = { 0, };
	equalhist[0] = hist[0];
	for (int i = 1; i < 256; i++) {
		equalhist[i] = equalhist[i - 1] + hist[i];
	}

	//정규화
	int norm_hist[256] = { 0, };
	for (int i = 0; i < 256; i++) {
		norm_hist[i] =  round(((double)equalhist[i] / (double)totalpixel) * 255);
	}
	
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			dst.at<uchar>(i, j) = norm_hist[src.at<uchar>(i, j) + 1];
		}
	}
	
	return dst;
}

Mat get_hist(Mat src)
{
	Mat dst = src;
	//히스토그램 구함
	int hist[256] = { 0, };
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			hist[dst.at<uchar>(i, j)] ++;
		}
	}
	return dst;
}


// spatial filter 적용
// k는 박스 크기. 
Mat Boxfilter(Mat src, int k)
{
	Mat dst(src.rows, src.cols, CV_8UC1, Scalar(0));
	int padder = (int)(k / 2);

	//제로패딩
	Mat zeropd = Mat::zeros(src.rows + padder + padder, src.cols + padder + padder, src.type());
	for (int row = 0; row < src.rows; row++) {
		for (int col = 0; col < src.cols; col++) {
			zeropd.at<uchar>(row + padder, col + padder) = src.at<uchar>(row, col);
		}
	}

	//convolution
	double dim = (double)1 / (double)(k * k);
	for (int row = padder; row < zeropd.rows - padder; row++) {
		for (int col = padder; col < zeropd.cols - padder; col++) {
			int sum = 0;
			for (int i = -padder; i <= padder; i++) {
				for (int j = -padder; j <= padder; j++) {
					sum += zeropd.at<uchar>(row + i, col + j);
				}
			}
			sum = round(dim * sum);
			dst.at<uchar>(row - padder, col - padder) = sum;
		}
	}

	//separable
	/*Mat dst2(src.rows, src.cols, CV_8UC1, Scalar(0));
	double sdim = (double)1 / (double)(k);
	for (int row = padder; row < zeropd.rows - padder; row++) {
		for (int col = padder; col < zeropd.cols - padder; col++) {
			int sum = 0;
			for (int i = -padder; i <= padder; i++) {
				sum += zeropd.at<uchar>(row + i, col);
			}
			sum = round(sdim * sum);
			zeropd.at<uchar>(row, col) = sum;
		}
	}
	for (int row = padder; row < zeropd.rows - padder; row++) {
		for (int col = padder; col < zeropd.cols - padder; col++) {
			int sum = 0;
			for (int i = -padder; i <= padder; i++) {
				sum += zeropd.at<uchar>(row, col + i);
			}
			sum = round(sdim * sum);
			dst2.at<uchar>(row - padder, col - padder) = sum;
		}
	}*/

	return dst;
}

// 가우시안의 분산은 각각 3, 7 (*6+1)
Mat Gaussianfilter(Mat src, int sigma)
{
	Mat dst(src.rows, src.cols, CV_8UC1, Scalar(0));
	int padder = (int)((sigma * 6 + 1) / 2); 
	Mat zeropd = Mat::zeros(src.rows + padder + padder, src.cols + padder + padder, src.type());

	// 가우시안 배열 생성
	double r, s = 2.0 * sigma * sigma;
	double sum = 0.0; // 정규화를 위한 총합
	double** kern = new double* [sigma * 6 + 1];
	for (int i = 0; i < sigma * 6 + 1; i++) {
		kern[i] = new double[sigma * 6 + 1];
	}
	for (int i = 0; i < sigma * 6 + 1; i++) {
		for (int j = 0; j < sigma * 6 + 1; j++) {
			kern[i][j] = 0;
		}
	}

	for (int x = -3 * sigma; x <= 3 * sigma; x++)
	{
		for (int y = -3 * sigma; y <= 3 * sigma; y++)
		{
			r = sqrt(x * x + y * y);
			kern[x + 3*sigma][y + 3*sigma] = (exp(-(r * r) / s)) / (M_PI * s);
			sum += kern[x + 3 * sigma][y + 3 * sigma];
		}
	}

	for (int i = 0; i < sigma * 6 + 1; i++) {
		for (int j = 0; j < sigma * 6 + 1; j++) {
			kern[i][j] /= sum;
		}
	}

	//임시 패딩된 배열 생성
	
	for (int row = 0; row < src.rows; row++) {
		for (int col = 0; col < src.cols; col++) {
			zeropd.at<uchar>(row + padder, col + padder) = src.at<uchar>(row, col);
		}
	}

	//convolution
	double dim = (double)1 / (double)((sigma * 6 + 1) * (sigma * 6 + 1));
	for (int row = padder; row < zeropd.rows - padder; row++) {
		for (int col = padder; col < zeropd.cols - padder; col++) {
			float ksum = 0;
			int yidx = 0, xidx = 0;
			for (int i = -padder; i <= padder; i++) {
				for (int j = -padder; j <= padder; j++) {
					ksum += zeropd.at<uchar>(row + i, col + j) * kern[xidx][yidx];
					xidx++;
				}
				xidx = 0;
				yidx++;
			}
			
			dst.at<uchar>(row - padder, col - padder) = round(ksum);
		}
	}

	delete kern;
	return dst;
}

Mat sobel(Mat src)
{
	//Mat dst = src;
	
	float sobeldata_x[] = { 1, 2, 1, 0,  0, 0, -1, -2, -1 };
	float sobeldata_y[] = { -1, 0, 1, -2, 0 ,2, -1, 0, 1 };

	Mat dst;
	Mat sobel_image1 = Mat::zeros(src.rows, src.cols, src.type());
	Mat sobel_image2 = Mat::zeros(src.rows, src.cols, src.type());

	//제로패딩
	Mat zeropd = Mat::zeros(src.rows + 2, src.cols + 2, src.type());

	for (int row = 0; row < src.rows; row++) {
		for (int col = 0; col < src.cols; col++) {
			zeropd.at<uchar>(row + 1, col + 1) = src.at<uchar>(row, col);
		}
	}
	
	//convolution
	for (int row = 1; row < zeropd.rows - 1; row++) {
		for (int col = 1; col < zeropd.cols - 1; col++) {
			int sum_x = 0;
			int idx = 0;
			for (int i = -1; i <= 1; i++) {
				for (int j = -1; j <= 1; j++) {
					sum_x += zeropd.at<uchar>(row + i, col + j) * sobeldata_x[idx];
					idx++;
				}
			}

			//0보다 작은 값은 0으로, 255보다 큰 값은 255로 채워줌
			if (sum_x < 0) { sum_x = 0; }
			else if (sum_x > 255) { sum_x = 255; }
			sobel_image1.at<uchar>(row - 1, col - 1) = (uchar)sum_x;
		}
	}

	for (int row = 1; row < zeropd.rows - 1; row++) {
		for (int col = 1; col < zeropd.cols - 1; col++) {
			int sum_y = 0;
			int idx = 0;
			for (int i = -1; i <= 1; i++) {
				for (int j = -1; j <= 1; j++) {
					sum_y += zeropd.at<uchar>(row + i, col + j) * sobeldata_y[idx];
					idx++;
				}
			}

			//0보다 작은 값은 0으로, 255보다 큰 값은 255로 채워줌
			if (sum_y <= 0) { sum_y = 0; }
			else if (sum_y >= 255) { sum_y = 255; }
			sobel_image2.at<uchar>(row - 1, col - 1) = (uchar)sum_y;
		}
	}

	add(sobel_image1, sobel_image2, dst);
	
	return dst;
}
Mat Avgpooling(Mat src, int k, int stride)
{
	int sz1 = ((src.rows) / stride);
	int sz2 = ((src.cols)/ stride);
	//cout << src.cols << ' ' << src.rows << ' ' <<
	//	sz2 << ' ' << sz1 << endl;
	int i = 0, j = 0;
	Mat dst(sz2, sz1, CV_8UC1, Scalar(0));
	for (i = 0; i < sz2 * k; i += stride) {
		for (j = 0; j < sz1 * k; j+=stride) {
			int kernelsum = 0;
			for (int x = 0; x < k; x++) {
				for (int y = 0; y < k; y++) {
					kernelsum += src.at<uchar>(i+x, j+y);
				}
			}
			kernelsum  /= (k * k);
			if (i >= sz2 * k || j >= sz1 * k)
			cout << i/stride << ' ' << j/stride << ' ' << kernelsum << endl;
			dst.at<uchar>(i / stride, j / stride) = round(kernelsum);
		}
	}

	return dst;
}

// https://preventionyun.tistory.com/32
Mat Downsampling(Mat src, int k1, int k2) {
	int sz1 = ((src.rows) / k1); //x
	int sz2 = ((src.cols) / k1); //y

	//cout << sz1 << " " << sz2 << endl;
	Mat ret(sz1, sz2, CV_8UC1, Scalar(256));

	double xscale = (double)sz1 / src.rows; 
	double yscale = (double)sz2 / src.cols; 

//	cout << scaleX << " " << scaleY << endl;
	for (int i = 0; i < src.cols; i++) { //y
		//int tmp1, tmp2;
		//double tmp3;
		for (int j = 0; j < src.rows; j++) { //x 
			int x = (int)(j * xscale); 
			int y = (int)(i * yscale); 
			ret.at<uchar>(x, y) = src.at<uchar>(j, i); 
			//tmp1 = y;
			//tmp2 = x;
		} 
		//cout << tmp1 << " "<<tmp2<<endl;
	} 
	
	return ret;
}

//https://m.blog.naver.com/PostView.naver?isHttpsRedirect=true&blogId=widebang&logNo=120210732324
// https://engineer-mole.tistory.com/14
// https://anothertechs.com/programming/cpp/opencv/opencv-low-high-pass-filter/

// https://m.blog.naver.com/PostView.naver?isHttpsRedirect=true&blogId=appcorn&logNo=220659361687

// 푸리에변환
Mat Freqfilter(Mat src, int flag) {
	if (flag == 0) {
		int optimalRows = getOptimalDFTSize(src.rows);
		int optimalCols = getOptimalDFTSize(src.cols);
		Mat padded;

		copyMakeBorder(src, padded, 0, optimalRows - src.rows, 0,
			optimalCols - src.cols, BORDER_CONSTANT, Scalar::all(0));

		Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
		Mat complexI;
		merge(planes, 2, complexI);         // Add to the expanded another plane with zeros
		dft(complexI, complexI);            // this way the result may fit in the source matrix

		// compute the magnitude and switch to logarithmic scale
		// => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
		split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
		magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude

		Mat magI = planes[0];
		magI += Scalar::all(1);                    // switch to logarithmic scale
		log(magI, magI);
		// crop the spectrum, if it has an odd number of rows or columns
		magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

		// rearrange the quadrants of Fourier image  so that the origin is at the image center
		int cx = magI.cols / 2;
		int cy = magI.rows / 2;

		Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
		Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
		Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
		Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
		Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
		q0.copyTo(tmp);
		q3.copyTo(q0);
		tmp.copyTo(q3);
		q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
		q2.copyTo(q1);
		tmp.copyTo(q2);
		normalize(magI, magI, 0, 1, NORM_MINMAX); // Transform the matrix with float values into a
												// viewable image form (float between values 0 and 1).

		return magI;
	}
	else {
		Mat outputImg;

		Mat planes[2] = { Mat_<float>(src.clone()), Mat::zeros(src.size(), CV_32F) };
		Mat complexI;
		merge(planes, 2, complexI);
		dft(complexI, complexI);
		split(complexI, planes);            // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
		planes[0].at<float>(0) = 0;
		planes[1].at<float>(0) = 0;
		// compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
		Mat imgPSD;
		magnitude(planes[0], planes[1], imgPSD);        //imgPSD = sqrt(Power spectrum density)
		pow(imgPSD, 2, imgPSD);                         //it needs ^2 in order to get PSD
		outputImg = imgPSD;

		return outputImg;
	}
}

Mat checker(Mat src, int cutoff, string str) {
	int sz1 = src.rows;
	int sz2 = src.cols;
	
	Mat ret = src;
	double ucenter = src.cols / 2, vcenter = src.rows / 2;

	if (str == "LPF") {
		for (int v = 0; v < src.rows; v++) {
			for (int u = 0; u < src.cols; u++) {
				double dist = sqrt(pow(u - ucenter, 2) + pow(v - vcenter, 2));

				Vec2f val = src.at<Vec2f>(v, u);
				if (dist <= cutoff) {
					val.val[0] *= 1.0;
					val.val[1] *= 1.0;
				}
				else {
					val.val[0] *= 0;
					val.val[1] *= 0;
				}
				ret.at<Vec2f>(v, u) = val;
			}
		}
	}
	else if (str == "Gauss") {
		for (int v = 0; v < src.rows; v++) {
			for (int u = 0; u < src.cols; u++) {
				double dist = sqrt(pow(u - ucenter, 2) + pow(v - vcenter, 2));

				// 2개짜리 배열 val
				Vec2f val = src.at<Vec2f>(v, u);
				val.val[0] *= exp(-dist * dist / (cutoff * cutoff * 2.0));
				val.val[1] *= exp(-dist * dist / (cutoff * cutoff * 2.0));

				ret.at<Vec2f>(v, u) = val;
			}
		}
	}

	//역변환
	Mat idftmat;
	dft(ret, idftmat, DFT_INVERSE | DFT_SCALE | DFT_REAL_OUTPUT);

	//위치재변경
	for (int y = 0; y < idftmat.rows; y++) {
		for (int x = 0; x < idftmat.cols; x++) {
			float fval = idftmat.at<float>(y, x);
			if ((x + y) % 2 == 1)
				idftmat.at<float>(y, x) = -fval;
		}
	}

	Mat dst;
	idftmat.convertTo(dst, CV_8U);
	//string title = str + to_string(cutoff);
	//imshow(title, dst);
	return dst;
}

//freq를 통해 들어온 푸리에 변환 도메인이 src. 이를 cutoff에 따라 분류해야 됨.
Mat LPF(Mat src, int cutoff) {
	Mat fimage = src;
	src.convertTo(fimage, CV_32F);

	// 위치 변경
	for (int y = 0; y < fimage.rows; y++) {
		for (int x = 0; x < fimage.cols; x++) {
			float fval = fimage.at<float>(y, x);
			if ((x + y) % 2 == 1)
				fimage.at<float>(y, x) = -fval;
		}
	}

	//시작. 일단 dft시킨다.
	Mat dftmat, dftmat_org;
	dft(fimage, dftmat_org, DFT_COMPLEX_OUTPUT);
	
	//실행
	dftmat_org.copyTo(dftmat);
	Mat ret = checker(dftmat, cutoff, "LPF");

	return ret;
}

Mat Gaussianfilterr(Mat src, int cutoff) {
	Mat fimage = src;
	src.convertTo(fimage, CV_32F);

	for (int y = 0; y < fimage.rows; y++) {
		for (int x = 0; x < fimage.cols; x++) {
			float fval = fimage.at<float>(y, x);
			if ((x + y) % 2 == 1)
				fimage.at<float>(y, x) = -fval;
		}
	}

	//시작
	Mat dftmat, dftmat_org;
	dft(fimage, dftmat_org, DFT_COMPLEX_OUTPUT);

	//실행.
	dftmat_org.copyTo(dftmat);
	//double D0[] = { 30.0,60.0,160.0,460.0 };
	//for (int i = 0; i < 1; i++) {
	Mat ret = checker(dftmat, cutoff, "Gauss");
	//}
	return ret;
}
/*
* 1. 일단 푸리에 필터를 구한다.
* 2. 해당 필터의 하얀 사각형 부분을 지워야됨.
* 2-1. 일단 파이썬 opencv로 추정 필수영역을 제외한 나머지를 다 검정색으로
* 2-2. 그 불필요한 부분을 찾은 다음, 그 부분의 일정 수 이상은 그냥 다 평균으로 메꾸기?
* 3. 그 지운 이미지를 푸리에 역변환.
*/

Mat notch(Mat imgIn) {
	imgIn.convertTo(imgIn, CV_32F);
	Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
	imgIn = imgIn(roi);

	// PSD calculation (start) - psd가 푸리에
	Mat imgPSD;
	imgPSD = Freqfilter(imgIn, 1);
	//calcPSD(imgIn, imgPSD);
	fftshift(imgPSD, imgPSD);
	normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);

	Mat H = Mat(roi.size(), CV_32F, Scalar(1));
	rectangle(H, Rect(Point(0, 0), Point(250, 250)), 0, -1);
	rectangle(H, Rect(Point(280, 0), Point(512, 250)), 0, -1);
	rectangle(H, Rect(Point(0, 280), Point(250, 512)), 0, -1);
	rectangle(H, Rect(Point(280, 280), Point(512, 512)), 0, -1);

	rectangle(H, Rect(Point(0, 250), Point(200, 280)), 0, -1);
	rectangle(H, Rect(Point(330, 250), Point(512, 280)), 0, -1);
	//H calculation (stop)
	
	// filtering (start)
	Mat imgOut;
	fftshift(H, H);
	filter2DFreq(imgIn, imgOut, H);
	// filtering (stop)

	imgOut.convertTo(imgOut, CV_8U);
	normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
	imwrite("result.tif", imgOut);
	imwrite("PSD.tif", imgPSD);
	fftshift(H, H);
	normalize(H, H, 0, 255, NORM_MINMAX);
	imwrite("filter.tif", H);

	return imgOut;
}


void fftshift(const Mat& inputImg, Mat& outputImg)
{
    outputImg = inputImg.clone();
    int cx = outputImg.cols / 2;
    int cy = outputImg.rows / 2;
    Mat q0(outputImg, Rect(0, 0, cx, cy));
    Mat q1(outputImg, Rect(cx, 0, cx, cy));
    Mat q2(outputImg, Rect(0, cy, cx, cy));
    Mat q3(outputImg, Rect(cx, cy, cx, cy));
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}

void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
{
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);
    dft(complexI, complexI, DFT_SCALE);
    Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
    Mat complexH;
    merge(planesH, 2, complexH);
    Mat complexIH;
    mulSpectrums(complexI, complexH, complexIH, 0);
    idft(complexIH, complexIH);
    split(complexIH, planes);
    outputImg = planes[0];
}

// Function calculates PSD(Power spectrum density) by fft with two flags
// flag = 0 means to return PSD
// flag = 1 means to return log(PSD)
//void calcPSD(const Mat& inputImg, Mat& outputImg, int flag)
//{
//    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
//    Mat complexI;
//    merge(planes, 2, complexI);
//    dft(complexI, complexI);
//    split(complexI, planes);            // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
//    planes[0].at<float>(0) = 0;
//    planes[1].at<float>(0) = 0;
//    // compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
//    Mat imgPSD;
//    magnitude(planes[0], planes[1], imgPSD);        //imgPSD = sqrt(Power spectrum density)
//    pow(imgPSD, 2, imgPSD);                         //it needs ^2 in order to get PSD
//    outputImg = imgPSD;
//}

이건 메인

#include "utils.h"

int main(int argc, char* argv[])
{
	/*
	//구현 1//

	Mat src_1 = imread("../images/pollen1.tif", IMREAD_ANYCOLOR);
	if (src_1.empty()) {
		cout << "Input image does not exist!" << endl;
		exit(-1);
	}
	int H = src_1.rows;
	int W = src_1.cols;

	Mat src_2 = imread("../images/pollen2.tif", IMREAD_ANYCOLOR);
	if (src_2.empty()) {
		cout << "Input image does not exist." << endl;
		exit(-1);
	}
	H = src_2.rows;
	W = src_2.cols;

	//imshow("Origin_1", src_1);
	//imshow("Origin_2", src_2);
	//a)
	Mat hist_1 = get_hist(src_1);
	Mat hist_2 = get_hist(src_2);

	//b) 
	Mat dst_1 = Histogram_equalization(src_1);
	//imshow("HE_1", dst_1);
	//imwrite("HE_1.tif", dst_1);

	Mat dst_2 = Histogram_equalization(src_2);
	//imshow("HE_2", dst_2);
	//imwrite("HE_2.tif", dst_2);


	//구현 2//

	src_2 = imread("../images/test_pattern.tif", IMREAD_ANYCOLOR);
	if (src_2.empty()) {
		cout << "Input image does not exist." << endl;
		exit(-1);
	}
	//a)

	Mat BF3x3 = Boxfilter(src_2, 3);
	Mat BF19x19 = Boxfilter(src_2, 19);
	Mat BF25x25 = Boxfilter(src_2, 25);
	//imshow("BF3x3", BF3x3);
	//imwrite("BF3x3.tif", BF3x3);

	//imshow("BF19x19", BF19x19);
	//imwrite("BF19x19.tif", BF19x19);

	//imshow("BF25x25", BF25x25);
	//imwrite("BF25x25.tif", BF25x25);

	//b)

	Mat GF3 = Gaussianfilter(src_2, 3);
	Mat GF7 = Gaussianfilter(src_2, 7);

	//imshow("GF3", GF3);
	//imwrite("GF3.tif", GF3);

	//imshow("GF7", GF7);
	//imwrite("GF7.tif", GF7);
	*/
	//구현3//

	Mat src_3 = imread("../images/contact-lens.tif", IMREAD_ANYCOLOR);
	if (src_3.empty()) {
		cout << "Input image does not exist." << endl;
		exit(-1);
	}
	
	//a)
	Mat AP = Avgpooling(src_3, 4, 4);
	//imshow("AP", AP);
	//imwrite("AP.tif", AP);

	//b)
	//Mat Sobel = sobel(src_3);
	//imshow("Sobel", Sobel);
	//imwrite("Sobel.tif", Sobel);

	//c)
	//Mat final = sobel(Avgpooling(src_3, 4, 4));
	//imshow("final", final);
	//imwrite("final.tif", final); 

	Mat barbara = imread("../images/barbara.tif", IMREAD_ANYCOLOR);
	Mat downsample = Downsampling(barbara, 3, 3);
	imshow("downsample", downsample);
	imwrite("downsample.tif", downsample);

	Mat freq = Freqfilter(barbara, 0); // 스펙트럼 출력
	imshow("Fourierbarbara", freq);
	imwrite("Fourierbarbara.tif", freq);

	int cutoff = 10;
	Mat tst = imread("../images/test_pattern.tif", IMREAD_ANYCOLOR);
	Mat tstforu = Freqfilter(tst, 0); // 스펙트럼 출력
	//imshow("Fouriertst", tstforu);
	imwrite("Fouriertst.tif", tstforu);

	Mat lpf = LPF(tst, cutoff); //ilpf 실행
	//Mat lpf = LPF(freq, cutoff);
	//imshow("LPF" + to_string(cutoff), lpf);
	imwrite("LPF" + to_string(cutoff) + ".tif", lpf);
	Mat gs = Gaussianfilterr(tst, cutoff); // 가우시안필터링
	//imshow("GAUSE" + to_string(cutoff), gs);
	imwrite("GAUSE" + to_string(cutoff) + ".tif", gs);
	
	Mat woman = imread("../images/noised_woman.tif", IMREAD_ANYCOLOR);
	Mat wmforu = Freqfilter(woman, 0); // 스펙트럼 출력
	imshow("Fourierwm", wmforu);
	imwrite("Fourierwm.tif", wmforu);
	
	Mat finale = notch(woman);
	imshow("finale", finale);
	imwrite("finale.tif", finale);

	//Mat w= imread("../images/noised_woman.jpg", IMREAD_GRAYSCALE);
	//Mat wm = Freqfilter(w); // 스펙트럼 출력
	//imshow("Fouriewm", wm);

	//Mat imgIn = imread("../images/period_input.jpg", IMREAD_GRAYSCALE);
	//Mat imgIn = imread("../images/noised_woman.tif", IMREAD_GRAYSCALE);
	//if (imgIn.empty()) //check whether the image is loaded or not
	//{
	//	cout << "ERROR : Image cannot be loaded..!!" << endl;
	//	return -1;
	//}
	//imgIn.convertTo(imgIn, CV_32F);
	//// it needs to process even image only
	//Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
	//cout << (imgIn.cols & -2) << " " << (imgIn.rows & -2) << endl; //640, 470
	//imgIn = imgIn(roi);
	//// PSD calculation (start) - psd가 푸리에
	//Mat imgPSD;
	//calcPSD(imgIn, imgPSD);
	//fftshift(imgPSD, imgPSD);
	//normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);
	//imwrite("PSD.tif", imgPSD);
	//imshow("psdd", imgPSD);
	//// PSD calculation (stop)

	////H calculation (start) --> 올바른 필터 좌표값 필요
	////rectangle( img , Rect , Scalar(b, g, r), thickness, lineType, shift )
	//Mat H = Mat(roi.size(), CV_32F, Scalar(1));
	//const int r = 21;
	//rectangle(H, Rect(Point(0, 0), Point(250, 250)), 0, -1);
	//rectangle(H, Rect(Point(280, 0), Point(512, 250)), 0, -1);
	//rectangle(H, Rect(Point(0, 280), Point(250, 512)), 0, -1);
	//rectangle(H, Rect(Point(280, 280), Point(512, 512)), 0, -1);

	//rectangle(H, Rect(Point(0, 250), Point(200, 280)), 0, -1);
	//rectangle(H, Rect(Point(330, 250), Point(512, 280)), 0, -1);

	//imshow("H-", H);
	////H calculation (stop)
	//
	//// filtering (start)
	//Mat imgOut;
	//fftshift(H, H);
	//filter2DFreq(imgIn, imgOut, H);
	//// filtering (stop)

	//imgOut.convertTo(imgOut, CV_8U);
	//normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
	//imshow("ret", imgOut);
	//imwrite("result.tif", imgOut);
	//fftshift(H, H);
	//normalize(H, H, 0, 255, NORM_MINMAX);
	//imwrite("filter.tif", H);

	//Mat fn = imread("C:/Users/is7se/source/repos/OpencvStudy/OpencvStudy/result.tif", IMREAD_ANYCOLOR);
	//Mat fnforu = Freqfilter(fn); // 스펙트럼 출력
	//imshow("abc", fnforu);


	waitKey(0);
	destroyAllWindows();

	return 0;
}