全变差TV去噪(split bregman) opencv

#include "subfunction.h" Mat psf2otf(Mat psf, Size size) { Mat otf = Mat::zeros(size.height, size.width, psf.type()); // Pad the PSF to outSize Size dftSize; // compute the size of DFT transform dftSize.width = getOptimalDFTSize(size.width); dftSize.height = getOptimalDFTSize(size.height); // allocate temporary buffers and initialize them with 0's Mat temp(dftSize, psf.type(), Scalar::all(0)); //copy psf to the top-left corners of temp Mat roipsf(temp, Rect(0, 0, psf.cols, psf.rows)); psf.copyTo(roipsf); // Circularly shift otf so that the "center" of the PSF is at the // (0,0) element of the array. Mat psf2 = temp.clone(); int cx = psf.cols / 2; int cy = psf.rows / 2; Mat p0(temp, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant Mat p1(temp, Rect(cx, 0, psf2.cols - cx, cy)); // Top-Right Mat p2(temp, Rect(0, cy, cx, psf2.rows - cy)); // Bottom-Left Mat p3(temp, Rect(cx, cy, psf2.cols - cx, psf2.rows - cy)); // Bottom-Right Mat q0(psf2, Rect(0, 0, psf2.cols - cx, psf2.rows - cy));// Top-Left - Create a ROI per quadrant Mat q1(psf2, Rect(psf2.cols - cx, 0, cx, psf2.rows - cy));// Top-Right Mat q2(psf2, Rect(0, psf2.rows - cy, psf2.cols - cx, cy)); // Bottom-Left Mat q3(psf2, Rect(psf2.cols - cx, psf2.rows - cy, cx, cy)); // Bottom-Right // swap quadrants (Top-Left with Bottom-Right) p0.copyTo(q3); p3.copyTo(q0); // swap quadrant (Top-Right with Bottom-Left) p1.copyTo(q2); p2.copyTo(q1); // Computer the OTF Mat planes[] = { Mat_<float>(psf2), Mat::zeros(psf2.size(), CV_32F) }; Mat complexI; merge(planes, 2, complexI); dft(complexI, complexI); otf = complexI; return otf(Range(0, size.height), Range(0, size.width)); } //max(abs(input) - 1 / lambda, 0).*sign(input) Mat softThreshold(Mat input, float lambda){ Mat output = Mat::zeros(input.size(), input.type()); Mat temp1 = Mat::zeros(input.size(), input.type()); Mat temp2 = Mat::zeros(input.size(), input.type()); temp1 = abs(input); temp1 = temp1 - (1 / lambda); temp1 = max(temp1, 0); temp2 = sign(input); output = temp1.mul(temp2); return output; } Mat sign(Mat input){ Mat output = Mat::zeros(input.size(), input.type()); Mat temp = abs(input); for (int i = 0; i < input.rows; i++){ for (int j = 0; j < input.cols; j++){ if (temp.at<float>(i, j) == 0) temp.at<float>(i, j) = 1; } } divide(input, temp, output); return output; } Mat fft2(Mat src) { int m = getOptimalDFTSize(src.rows); int n = getOptimalDFTSize(src.cols); Mat padded; copyMakeBorder(src, padded, 0, m - src.rows, 0, n - src.cols, BORDER_CONSTANT, Scalar::all(0)); Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32FC1) }; Mat complexI; merge(planes, 2, complexI); dft(complexI, complexI); return complexI; } Mat diffX(Mat input){ Mat Ux = Mat::zeros(input.size(), input.type()); for (int i = 0; i < input.cols; i++) { if (i == input.cols - 1) Ux.col(i) = input.col(0) - input.col(i); else Ux.col(i) = input.col(i + 1) - input.col(i); } return Ux; } Mat diffY(Mat input){ Mat Uy = Mat::zeros(input.size(), input.type()); for (int i = 0; i < input.rows; i++) { if (i == input.cols - 1) Uy.row(i) = input.row(0) - input.row(i); else Uy.row(i) = input.row(i + 1) - input.row(i); } return Uy; } Mat conj(Mat input){ Mat output; Mat phane[] = { Mat::zeros(input.size(), input.type()), Mat::zeros(input.size(), input.type()) }; split(input, phane); phane[1] = 0 - phane[1]; merge(phane, 2, output); return output; } Mat complexMul(Mat src1, Mat src2){ Mat temp[] = { Mat::zeros(src1.size(), src1.type()), Mat::zeros(src1.size(), src1.type()) };; Mat temp1[] = { Mat::zeros(src1.size(), src1.type()), Mat::zeros(src1.size(), src1.type()) }; Mat temp2[] = { Mat::zeros(src1.size(), src1.type()), Mat::zeros(src1.size(), src1.type()) }; Mat temp3, temp4, temp5, temp6, output; split(src1, temp1); split(src2, temp2); multiply(temp1[0], temp2[0], temp3); multiply(temp1[1], temp2[1], temp4); temp[0] = temp3 - temp4; multiply(temp1[0], temp2[1], temp5); multiply(temp1[1], temp2[0], temp6); temp[1] = temp5 + temp6; merge(temp, 2, output); return output; } void displayImage(Mat Image, string winName){ Image = norm_0_255(Image); namedWindow(winName, CV_WINDOW_AUTOSIZE); imshow(winName, Image); } Mat norm_0_255(Mat src) { Mat dst; switch (src.channels()){ case 1: normalize(src, dst, 0, 255, CV_MINMAX, CV_8UC1); break; case 3: normalize(src, dst, 0, 255, CV_MINMAX, CV_8UC3); break; default: src.copyTo(dst); break; } return dst; } double PSNR(Mat &src1, Mat &src2){ Mat temp1, temp2; double a = 0; if (src1.isContinuous()) src1.reshape(1, src1.total()).copyTo(temp1); else{ src1.copyTo(temp1); temp1.reshape(1, src1.total()); } if (src2.isContinuous()) src2.reshape(1, src2.total()).copyTo(temp2); else{ src1.copyTo(temp2); temp2.reshape(1, src2.total()); } temp1 = temp1 - temp2; pow(temp1, 2, temp1); temp1 = mean(temp1); sqrt(temp1, temp1); a = temp1.at<float>(0, 0); a = 255 / a; a = 20 * log10(a); return a; }