nlm算法matlab转c程序

admin_tmp

已于 2023-08-02 09:35:38 修改

阅读量124

点赞数 1

文章标签： c语言

于 2023-08-02 00:24:26 首次发布

本文链接：https://blog.csdn.net/weixin_44745808/article/details/132051696

版权

nlm-C程序

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"

void mirr_padding_h(unsigned char* src, unsigned char* dst, int src_w, int src_h, int pad_num)
{
	int dst_w = src_w;
	int dst_h = src_h+(2*pad_num);
	int *idx1;
	int *idx2;
	idx1 = malloc(pad_num * sizeof(int));
	idx2 = malloc(pad_num * sizeof(int));
	for (int i = 0; i < pad_num; i++)
		idx1[i] = pad_num - 1 - i;
	for (int h = 0; h<pad_num; h++) {
		for(int w=0; w<dst_w; w++){
			dst[w + dst_w * h] = src[w + src_w * idx1[h]];
		}
	}
	for(int h=pad_num; h<src_h+pad_num; h++){
		for(int w=0; w<dst_w; w++){
			dst[w + dst_w * h] = src[w + src_w * (h-pad_num)];
		}
	}
	for (int i = 0; i < pad_num; i++) {
		idx2[i] = src_h - 1 - i;
		
	}
	for(int h=src_h+pad_num; h<dst_h; h++){
		for(int w=0; w<dst_w; w++){
			dst[w + dst_w * h] = src[w + src_w * idx2[h - src_h - pad_num]];
		}
	}
	free(idx1);
	free(idx2);
}
void mirr_padding_w(unsigned char* src, unsigned char* dst, int src_w, int src_h, int pad_num)
{
	int dst_w = src_w+(2*pad_num);
	int dst_h = src_h;
	int *idx1;
	int *idx2;
	idx1 = malloc(pad_num * sizeof(int));
	idx2 = malloc(pad_num * sizeof(int));
	for (int i = 0; i < pad_num; i++){
		idx1[i] = pad_num - 1 - i;
		idx2[i] = src_w - 1 - i;
		
	}
	for (int h = 0; h<dst_h; h++) {
		for(int w=0; w<pad_num; w++){
			dst[w + dst_w * h] = src[idx1[w] + src_w * h];
		}
		for(int w=pad_num; w<src_w+pad_num; w++){
			dst[w + dst_w * h] = src[w-pad_num + src_w * h];
		}
		for(int w=src_w+pad_num; w<dst_w; w++){
			dst[w + dst_w * h] = src[idx2[w-src_w-pad_num] + src_w * h];
		}
	}
	free(idx1);
	free(idx2);
}
int padding(unsigned char* src, unsigned char* dst, int padding_width, int padding_high, int width, int high, int ds, int Ds)
{
	int idx_mirror;
	for (int h = 0; h < ds + Ds; h++) {
		idx_mirror = (Ds + ds - h) * width;
		for (int w = 0; w < width; w++){
			dst[w + h * padding_width + ds + Ds] = src[w + idx_mirror];
		}
	}
	for (int h = Ds + ds + high; h < padding_high; h++) {
		idx_mirror = width * (high - 1 - h + Ds + ds + high);
		for (int w = 0; w < width; w++){
			dst[w + h * padding_width + ds + Ds] = src[w + idx_mirror];
		}
	}

	for (int h = 0; h < high; h++){
		int idx = h * width;
		for (int w = 0; w < width; w++){
			dst[(ds + Ds + h) * padding_width + ds + Ds + w] = src[w + idx];
		}
	}

	for (int h = 0; h < padding_high; h++){
		for (int w = 0; w < ds + Ds; w++){
			dst[h * padding_width + w] = dst[h * padding_width + (ds + Ds - w - 1) + ds + Ds];
		}
	}
	for (int h = 0; h < padding_high; h++){
		for (int w = width + Ds + ds; w < padding_width; w++){
			dst[h * padding_width + w] = dst[h * padding_width + ds + Ds + width - 1 - w + width + Ds + ds];
		}
	}
	return 0;
}
unsigned char* mat_img(unsigned char* src, int src_w, int st_w, int sp_w, int st_h, int sp_h)
{
	int dst_w = sp_w - st_w + 1;
	int dst_h = sp_h - st_h + 1;
	unsigned char* dst = NULL;
	dst = (unsigned char*)malloc(dst_w * dst_h * sizeof(unsigned char));
	for (int h = 0; h < dst_h; h++) {
		for (int w = 0; w < dst_w; w++) {
			int tmp = src[w + st_w + (h + st_h) * src_w];
			dst[w + h * dst_w] = src[w + st_w + (h + st_h) * src_w];
		}
	}
	return dst;
}
int read_file(unsigned char* data, int width, int high)
{
	FILE* fp = NULL;
	int count;
	fp = fopen("E:/Projects/matlab/src.y", "rb");
	if (fp == NULL)
	{
		printf("read file failed\n");
		return 2;
	}
	count = fread(data, sizeof(unsigned char), width * high, fp);
	fclose(fp);
	return 0;
}
int main()
{
	int Ds = 5;
	int ds = 2;
	int h_w = 10;
	int img_w = 128;
	int img_h = 128;

	int simg_w = img_w + ds;
	int simg_h = img_h + ds;
	int pad_w = img_w + (ds + Ds) * 2;
	int pad_h = img_h + (ds + Ds) * 2;

	unsigned char* src;
	unsigned char* src_padh;
	unsigned char* PaddedImg;
	unsigned char* image;
	unsigned char* wimage;
	int* diff;
	int* diff_sum1;
	int* diff_sum2;
	double* distance;
	double* weight;
	double* sumimage;
	double* sumweight;
	double* maxweight;
	double* dst;
	src = (unsigned char*)malloc(img_w * img_h * sizeof(unsigned char));
	src_padh = (unsigned char*)malloc(img_w * (img_h + (Ds + ds)*2) * sizeof(unsigned char));
	PaddedImg = (unsigned char*)malloc(pad_w * pad_h * sizeof(unsigned char));
	image = (unsigned char*)malloc(simg_w * simg_h * sizeof(unsigned char));
	wimage = (unsigned char*)malloc(simg_w * simg_h * sizeof(unsigned char));
	diff = (int*)malloc(simg_w * simg_h * sizeof(int));
	diff_sum1 = (int*)malloc(simg_w * simg_h * sizeof(int));
	diff_sum2 = (int*)malloc(simg_w * simg_h * sizeof(int));
	distance = (double*)malloc(img_w * img_h * sizeof(double));
	weight = (double*)malloc(img_w * img_h * sizeof(double));
	sumimage = (double*)malloc(img_w * img_h * sizeof(double));
	sumweight = (double*)malloc(img_w * img_h * sizeof(double));
	maxweight = (double*)malloc(img_w * img_h * sizeof(double));
	dst = (double*)malloc(img_w * img_h * sizeof(double));
	read_file(src, img_w, img_h);
	//padding(src,PaddedImg,pad_w,pad_h,img_w,img_h,ds,Ds);
	mirr_padding_h(src, src_padh, img_w, img_h, ds + Ds);
	mirr_padding_w(src_padh, PaddedImg, img_w, img_h+(ds+Ds)*2, ds + Ds);
	int debug = 1;
	image = mat_img(PaddedImg, pad_w, Ds, Ds + img_w + ds - 1, Ds, Ds + img_h + ds - 1);
	for (int r = -Ds; r <= Ds; r++) {
		for (int s = -Ds; s <= Ds; s++) {
			// 跳过当前点偏移
			if (r == 0 && s == 0)
				continue;
			
			// 求得差值积分图
			for (int h = 0; h < simg_h; h++) {
				for (int w = 0; w < simg_w; w++) {
					wimage[w + h * simg_w] = PaddedImg[w + Ds + s + (h + Ds + r) * pad_w];
				}
			}
			for (int h = 0; h < simg_h; h++) {
				for (int w = 0; w < simg_w; w++) {
					diff[w + h * simg_w] = ((int)image[w + h * simg_w] - (int)wimage[w + h * simg_w]) * ((int)image[w + h * simg_w] - (int)wimage[w + h * simg_w]);
				}
			}
			// 积分图
			for (int h = 0; h < simg_h; h++) {
				for (int w = 0; w < simg_w; w++) {
					if (w == 0)
						diff_sum1[h * simg_w] = diff[h * simg_w];
					else
						diff_sum1[w + h * simg_w] = diff_sum1[(w - 1) + h * simg_w] + diff[w + h * simg_w];
				}
			}

			for (int w = 0; w < simg_w; w++) {
				for (int h = 0; h < simg_h; h++) {
					if (h == 0)
						diff_sum2[w] = diff_sum1[w];
					else
						diff_sum2[w + h * simg_w] = diff_sum2[w + (h - 1) * simg_w] + diff_sum1[w + h * simg_w];
				}
			}

			// 计算距离
			int tmp;
			for (int h = 0; h < img_h; h++) {
				for (int w = 0; w < img_w; w++) {
					tmp = diff_sum2[simg_w - img_w + w + (h + simg_h - img_h) * simg_w]
						+ diff_sum2[w + h * simg_w]
						- diff_sum2[simg_w - img_w + w + h * simg_w]
						- diff_sum2[w + (h + simg_h - img_h) * simg_w];
					distance[w + h * img_w] = (double)tmp / (double)((2 * ds + 1) * (2 * ds + 1));
					weight[w + h * img_w] = exp(-distance[w + h * img_w] / (double)(h_w * h_w));
				}
			}
			
			// 计算权重并获得单个偏移下的加权图像
			for (int h = 0; h < img_h; h++) {
				for (int w = 0; w < img_w; w++) {
					if (r == -Ds && s == -Ds) {
						sumimage[w + h * img_w] = weight[w + h * img_w] * (double)wimage[w + ds + (h + ds) * simg_w];
						sumweight[w + h * img_w] = weight[w + h * img_w];
						maxweight[w + h * img_w] = weight[w + h * img_w];
					}
					else {
						sumimage[w + h * img_w] = sumimage[w + h * img_w] + weight[w + h * img_w] * (double)wimage[w + ds + (h + ds) * simg_w];
						sumweight[w + h * img_w] = sumweight[w + h * img_w] + weight[w + h * img_w];
						maxweight[w + h * img_w] = (maxweight[w + h * img_w] > weight[w + h * img_w]) ? maxweight[w + h * img_w] : weight[w + h * img_w];
					}
				}
			}
			/*if (s == Ds-1) {
				for (int i = 0; i < 1440; i++) {
					printf("%4d :%f\n", i + 1, sumimage[i * 1440]);
				}
				int debug = 1;
			}*/
		}

	}


	for (int h = 0; h < img_h; h++) {
		for (int w = 0; w < img_w; w++) {
			sumimage[w + h * img_w] = sumimage[w + h * img_w] + maxweight[w + h * img_w] * image[w + ds + (h + ds) * simg_w];
			sumweight[w + h * img_w] = sumweight[w + h * img_w] + maxweight[w + h * img_w];
			dst[w + h * img_w] = sumimage[w + h * img_w] / sumweight[w + h * img_w];
		}
	}


	// 滤波结果输出
	FILE* fp = NULL;
	size_t count = 0;
	fp = fopen("E:/Projects/matlab/yuv420_nlm.y", "wb");
	if (fp != NULL) {
		count = fwrite(dst, sizeof(double), (size_t)img_w * (size_t)img_h, fp);
	}
	else
		return -1;
	if (count != (size_t)img_w * (size_t)img_w) {
		printf("write failed\n");
	}
	fclose(fp);

	free(src);
	free(src_padh);
	free(PaddedImg);
	free(image);
	free(wimage);
	free(diff);
	free(diff_sum1);
	free(diff_sum2);
	free(distance);
	free(weight);
	free(sumimage);
	free(sumweight);
	free(maxweight);
	free(dst);
	return 0;
}

MATALB检查C计算结果程序

clc
w=128;
h=128;
fp = fopen("yuv420_nlm.y",'r');
I= fread(fp,w*h,"double");
fclose(fp);
demo=zeros(h,w);
for i=1:h
    for j=1:w
        demo(i,j) = I(j+(i-1)*w);
    end
end

err = demo-dst_1;
if max(max(err))<0.01
    fprintf(' Success!\n Max error is %d\n',max(max(err)))
else
    fprintf(' Failed!\n Max error is %d\n',max(max(err)))
end

MATLAB nlm算法 + 导出素材

clc
clear all
close all
 
%---------------------------%
%       input
%---------------------------%
src=imread('E:\Projects\matlab\nlm\Non_Local_Means_Denoise-main\img\noise\gaussian.png');
src=rgb2gray(src);
src=double(src);
figure,imshow(src,[]),title('src image')

%---------------------------%
%       写入文件
%---------------------------%
fp = fopen("src.y",'w');
fwrite(fp,src',"uint8");
fclose(fp); 

%---------------------------%
%       parameter
%---------------------------%
[m,n]=size(src);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(src,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像，便于处理
 
%---------------------------%
%  Non-Local Means Denoising
%---------------------------%
sumimage=zeros(m,n);
sumweight=zeros(m,n);
maxweight=zeros(m,n);
image=PaddedImg(1+Ds:Ds+m+ds,1+Ds:Ds+n+ds);
[M,N]=size(image);
for r=-Ds:Ds
    for s=-Ds:Ds
        
        %跳过当前点偏移
        if(r==0&&s==0)
            continue;
        end
        
        %求得差值积分图
        wimage=PaddedImg(1+Ds+r:Ds+m+ds+r,1+Ds+s:Ds+n+ds+s);
%         wimage=image(r:M+r,s:N+s);
        diff=image-wimage;
        diff=diff.^2;
        J=cumsum(diff,1);
        J=cumsum(J,2);
        
        %计算距离
        distance=J(M-m+1:M,N-n+1:N)+J(1:m,1:n)-J(M-m+1:M,1:n)-J(1:m,N-n+1:N);
        distance=distance/((2*ds+1).^2);
        
        %计算权重并获得单个偏移下的加权图像
        weight=exp(-distance./(h*h));
        sumimage=sumimage+weight.*wimage(ds+1:ds+m,ds+1:ds+n);
        sumweight=sumweight+weight;
        maxweight=max(maxweight,weight);
        if s==Ds
            debug = 1;
        end
    end
    if r==-Ds+1
        debug = 1;
    end
end
sumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n);
sumweight=sumweight+maxweight;
dst_1=sumimage./sumweight;
 
%---------------------------%
%       output
%---------------------------%
figure,imshow(dst_1,[]),title('dst');