图像处理笔记

阅读数:0次

将plot后的图像保存到本地(matlab)

参考
https://blog.csdn.net/qq_16775293/article/details/95319054

黑白图像灰度扫描

要求：提取某一行或者某一列的像素矢量，并呈现扫描结果
process.m 代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
foldername = 'assignment01_images\';
filename = 'cameraman.tif';
[a,b,c] = fileparts([foldername,filename]);
img = imread([foldername,filename]);
[dm, dn] = size(img);
sr = scanLine4e(img, fix((dm+1)/2), 'row');
sc = scanLine4e(img, fix((dn+1)/2), 'column'); 
h=figure(1);
plot(sr, 'r')
title([b, '图像灰度扫描']);
legend('中心行的像素矢量');
hold on
h2 = figure(2);
plot(sc, 'b')
title([b, '图像灰度扫描']);
legend('中心列的像素矢量');
print(h,'-djpeg',[b, '_row.jpg']);
print(h2,'-djpeg',[b, '_col.jpg']);

scanLine4e.m 代码

1
2
3
4
5
6
7
8
9
10
function s = scanLine4e(f, I, loc)
if(nargin<3)
   loc = 'row';
end
if strcmp(loc,'row')
    s = f(I, :);
end 
if strcmp(loc,'column')
    s = f(:, I);
end

彩色图像转换为黑白图

  1. 基础

    • 最大值法：使R,G,B的值等于3值中最大的一个，R=G=B=max（R,G,B），最大值法会形成亮度很高的灰度图像
    • 平均最突出和最不突出的部分：Y=(max(R,G,B)+min(R,G,B))/2
    • 平均值法：Y=(R+G+B)/3，会形成较为柔和的灰度图像
    • 加权平均法： 人对绿色最敏感，Y=0.2989∗R+0.5870∗G+0.1140∗B，会得到最合理的灰度图像

  2. 代码
    process.m 代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
     foldername = 'assignment01_images\';
     filename = 'mandril_color.tif';
     [a,b,c] = fileparts([foldername,filename]);
     img = imread([foldername,filename]);
     %R=img;R(:,:,[2 3])=0;
     %G=img;G(:,:,[1 3])=0;
     %B=img;B(:,:,[1 2])=0;
     %figure,imshow(R,[]), title('R分量'); 
     %figure,imshow(G,[] ), title('G分量'); 
     %figure,imshow(B,[]), title('B分量'); 
     g1 = rgb1gray(img, 'average');
     figure,imshow(g1,[]), title([b,' average'],'Interpreter','none');
     g2 = rgb1gray(img, 'NTSC');
     figure,imshow(g2,[]), title([b,' NTSC'],'Interpreter','none');
     imwrite(g1,strcat(b,'_average.jpg'));
     imwrite(g2,strcat(b,'_NTSC.jpg'));
    ```  
    rgb1gray.m 代码
    ``` matlab
     function g = rgb1gray(f, method)
     if(nargin<2)
         method = 'NTSC';
     end
     R=f(:,:,1);
     G=f(:,:,2);
     B=f(:,:,3);
     if strcmp(method,'average')
         g = R/3+G/3+B/3;
     end 
     if strcmp(method, 'NTSC')
         g = R*0.2989+ G*0.5870+ B*0.1140;
     end

  3. 参考
    https://blog.csdn.net/weixin_41645983/article/details/90641408
    https://blog.csdn.net/giantpoplar/article/details/47295979
    https://blog.csdn.net/caomin1hao/article/details/81876683
    https://blog.csdn.net/lihe4151021/article/details/89366425
    http://www.luyixian.cn/news_show_3756.aspx

二维卷积实现

  1. 原理
    主要步骤 卷积核旋转180度->填充padding->卷积
    卷积核旋转180度原因：
    https://my.oschina.net/datadev/blog/1859197
    https://blog.csdn.net/tina_ttl/article/details/51034809
    注意！！matlab下标从1开始

  2. 坑点
    注意数据类型和精度
    当卷积核输入double类型精度不同时，twoConv conv2 imfilter 即使同样是double类型，输出的精度也不一定相同。这里误差取 <1e-8 为相同
    imfilter 输出和输入类型一样，要自己显式转double，conv2 自动转double

  3. 代码
    卷积函数实现：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    function g = twodConv(f, w, padding)
    if(nargin < 3)
      padding = 'zero';
     end
     f = double(f);
     w = double(w);
     [fm, fn] = size(f);
     [wm, wn] = size(w);
     if strcmp(padding, 'replicate')
         t = [f(:,1)*ones(1, (wn-1)/2) f f(:,end)*ones(1,(wn-1)/2)];
         t = [ones((wm-1)/2, 1)*t(1,:); t; ones((wm-1)/2, 1)*t(end, :)];
     end
     if strcmp(padding, 'zero')
         t = zeros(fm+wm-1, fn+wn-1);
         t(((wm-1)/2+1):((wm-1)/2+fm), ((wn-1)/2+1):((wn-1)/2+fn)) = f;
     end
     g = zeros(fm,fn);
     w = rot90(w,2);
     %----method 1
     for i = 1: fm
         for j = 1: fn
             g(i, j) = sum(sum(t(i:i+wm-1, j:j+wn-1).*w));
         end
     end
     %----method 2
     %for i = 1: fm
     %    for j = 1: fn
     %        for x = 1:wm
     %            for y = 1:wn
     %                g(i,j) = g(i,j) + t(i+x-1, j+y-1)* w(x,y);
     %            end
     %        end
     %    end
     %end
     %----method 3
     %for i = ((wm-1)/2+1):((wm-1)/2+m)
     %    for j = ((wn-1)/2+1):((wn-1)/2+n)
     %       for x = 1:wm
     %           for y = 1:wn
     %               %[i-(wm-1)/2, j-(wn-1)/2, i-(wm-1)/2-1+x, j-(wn-1)/2-1+y, x, y]
     %               g(i-(wm-1)/2,j-(wn-1)/2) =  g(i-(wm-1)/2,j-(wn-1)/2) + t(i-(wm-1)/2-1+x, j-(wn-1)/2-1+y)*w(x,y);
     %           end
     %       end
     %    end
     %end

    主函数实现:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    %-----read data
     foldername = 'images\';
     %filename = 'einstein.tif';
     [a,b,c] = fileparts([foldername,filename]);
     img = imread([foldername,filename]);
     %img = im2double(img);

     %-----choose kernel
     kernel = [0,1,0; 1,-4,1; 0,1,0];
     %kernel = randi([-10,10], 3);
     %kernel = 10*rand(3);

     %-----call method
     z0 = twodConv(img, kernel);
     z1 = conv2(double(img), double(kernel),'same');
     z2 = imfilter(double(img), double(kernel), 'conv', 0, 'same');
     r0 = twodConv(img, kernel, 'replicate');
     r1 = imfilter(double(img), double(kernel), 'conv', 'replicate', 'same');

     %----show
     h1=figure;imshow(img,[]), title([b,'_origin'],'Interpreter','none');
     h2=figure;imshow(z0,[]), title([b,'_twodConv_zero'],'Interpreter','none');
     h3=figure;imshow(z1,[]), title([b,'_conv2_zero'],'Interpreter','none');
     h4=figure;imshow(z2,[]), title([b,'_imfilter_zero'],'Interpreter','none');
     h5=figure;imshow(r0,[]), title([b,'_twodConv_replicate'],'Interpreter','none');
     h6=figure;imshow(r1,[]), title([b,'_imfilter_replicate'],'Interpreter','none');

     %----compare
     % diff:num, num不为0说明存在不同，num表示不同的像素个数
     EPS = 1e-9;
     fprintf('z0=z1? diff:%d\n',sum(sum(abs(z1-z0)>EPS)))
     fprintf('z0=z2? diff:%d\n',sum(sum(abs(z2-z0)>EPS)))
     fprintf('z1=z2? diff:%d\n',sum(sum(abs(z1-z2)>EPS)))
     fprintf('r0=r1? diff:%d\n',sum(sum(abs(r1-r0)>EPS)))
     fprintf('z2=r1? diff:%d\n',sum(sum(abs(r1-z2)>EPS)))

     %----compare my vs sys
     Ip=uint8(z1)-uint8(z0);
     h7=figure;imshow(Ip,[]), title([b,'_my_sys_uint8_diff'],'Interpreter','none');
     Ip2 = imcomplement(Ip);
     h8=figure;imshow(Ip2,[]), title([b,'_my_sys_uint8_flip_diff'],'Interpreter','none');
     Ip3=z1-z0;
     h9=figure;imshow(Ip3,[]), title([b,'_my_sys_double_diff'],'Interpreter','none');
     Ip4 = imcomplement(Ip3);
     h10=figure;imshow(Ip4,[]), title([b,'_my_sys_double_flip_diff'],'Interpreter','none');

     %----compare zero vs replicate
     Ip5=r0-z0;
     %Ip=uint8(r0)-uint8(z0);
     h11=figure;imshow(Ip5,[]), title([b,'_diff'],'Interpreter','none');
     Ip6 = imcomplement(Ip5);
     h12=figure;imshow(Ip6,[]), title([b,'_flip_diff'],'Interpreter','none');

     %----save
     resultsFolder='results\question1\';
     print(h2,'-djpeg',[resultsFolder, b, '_zero_twoConv.jpg']);
     print(h3,'-djpeg',[resultsFolder, b, '_zero_conv2.jpg']);
     print(h4,'-djpeg',[resultsFolder, b, '_zero_imfilter.jpg']);
     print(h5,'-djpeg',[resultsFolder, b, '_replicate_twoConv.jpg']);
     print(h6,'-djpeg',[resultsFolder, b, '_replicate_imfilter.jpg']);
     print(h7,'-djpeg',[resultsFolder, b, '_my_sys_uint8_diff.jpg']);
     print(h8,'-djpeg',[resultsFolder, b, '_my_sys_uint8_flip_diff.jpg']);
     print(h9,'-djpeg',[resultsFolder, b, '_my_sys_double_diff.jpg']);
     print(h10,'-djpeg',[resultsFolder, b, '_my_sys_double_flip_diff.jpg']);
     print(h11,'-djpeg',[resultsFolder, b, '_diff.jpg']);
     print(h12,'-djpeg',[resultsFolder, b, '_flip_diff.jpg']);

  4. 验证

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
       fprintf('img2: %f %f\n',min(img2(:)), max(img2(:)) 
       imequal   %整数类型
       sum(sum(abs(r1-r0)>EPS) %浮点类型
       % 不同像素的定位
       [row,col,dim] = ind2sub(size(img3), find(img3 ~= img31)); 
       locations2 = [row, col, dim];
       disp(locations2);
       ```  
    5. 参考  
       matlab文档：  
       https://ww2.mathworks.cn/help/matlab/ref/conv2.html#bvhtlhr-1  
       图像增强一些知识：   
       https://blog.csdn.net/u012679707/article/details/78220298  
       边界填充(replicate代码实现)：   
       https://blog.csdn.net/lanchunhui/article/details/51419471  
       卷积代码：  
       https://blog.csdn.net/qq_38327353/article/details/80145276  
       imfilter使用：  
       https://blog.csdn.net/zhuwei0710/article/details/68169317  
       conv2 imfilter区别：  
       https://blog.csdn.net/carson2005/article/details/49640885  
       https://blog.csdn.net/qq_35860352/article/details/84025936   
       常用算子：  
       https://blog.csdn.net/Chaolei3/article/details/79809703 

    ## 归一化二维高斯滤波
    1. 代码  
       高斯核代码实现：
       ``` matlab
       function w = gaussKernel(sig, m)
        if(nargin < 2)
            m = ceil(3*sig)*2 + 1;  
        end
        if m < ceil(3*sig)*2 + 1
            fprintf('m is too small, please input again!')
            return
        else
            w = zeros(m, m);
            center = fix(m/2)+1; %注意取整
            for i = 1:m
                for j = 1: m
                    temp = double((i-center)^2+(j-center)^2);
                    w(i, j) = exp(-temp/(2*sig*sig))/(2*pi*sig*sig);
                end
            end
            w = w./sum(w(:));
        end

    主函数实现：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    %-----read data
     foldername = 'images\';
     filename = 'einstein.tif';
     [a,b,c] = fileparts([foldername,filename]);
     img = imread([foldername,filename]);
     %img = im2double(img);

     %-----kernel
     sig = 5;
     m = ceil(3*sig)*2 + 1;
     sysKernel = fspecial('gaussian',[m, m],sig);
     myKernel = gaussKernel(sig, m);
     EPS = 1e-9;
     fprintf('sysKernel=myKernel? diff:%d\n',sum(sum(abs(sysKernel-myKernel)>EPS)));

     %-----call method
     z0 = twodConv(double(img),double(myKernel));
     z1 = imfilter(double(img),double(myKernel),0,'same');
     z2 = imfilter(double(img),double(sysKernel),0,'same');
     z3 = imgaussfilt(double(img), sig,'FilterSize',m,'padding',0);
     r0 = twodConv(double(img),double(myKernel), 'replicate');
     r1 = imfilter(double(img),double(myKernel),'replicate','same');
     r2 = imfilter(double(img),double(sysKernel),'replicate','same');
     r3 = imgaussfilt(double(img), sig,'FilterSize',m,'padding','replicate');

     %----show
     h0=figure;imshow(img,[]), title([b,'_origin'],'Interpreter','none');
     h1=figure;imshow(z0,[]), title([b,'_my_twodConv_zero_sig_', num2str(sig)],'Interpreter','none');
     h2=figure;imshow(z1,[]), title([b,'_my_imfilter_zero_sig_', num2str(sig)],'Interpreter','none');
     h3=figure;imshow(z2,[]), title([b,'_sys_imfilter_zero_sig_', num2str(sig)],'Interpreter','none');
     h4=figure;imshow(z3,[]), title([b,'_sys_imgaussfilt_zero_sig_', num2str(sig)],'Interpreter','none');
     h5=figure;imshow(r0,[]), title([b,'_my_twodConv_replicate_sig_', num2str(sig)],'Interpreter','none');
     h6=figure;imshow(r1,[]), title([b,'_my_imfilter_replicate_sig_', num2str(sig)],'Interpreter','none');
     h7=figure;imshow(r2,[]), title([b,'_sys_imfilter_replicate_sig_', num2str(sig)],'Interpreter','none');
     h8=figure;imshow(r3,[]), title([b,'_sys_imgaussfilt_replicate_sig_', num2str(sig)],'Interpreter','none');

     %----compare
     % diff:num, num不为0说明存在不同，num表示不同的像素个数
     EPS = 1e-9;
     fprintf('z0=z1? diff:%d\n',sum(sum(abs(z0-z1)>EPS)))
     fprintf('z0=z2? diff:%d\n',sum(sum(abs(z0-z1)>EPS)))
     fprintf('z0=z3? diff:%d\n',sum(sum(abs(z0-z3)>EPS)))
     fprintf('z1=z2? diff:%d\n',sum(sum(abs(z1-z2)>EPS)))
     fprintf('z1=z3? diff:%d\n',sum(sum(abs(z1-z3)>EPS)))
     fprintf('z2=z3? diff:%d\n',sum(sum(abs(z2-z3)>EPS)))

     fprintf('r0=r1? diff:%d\n',sum(sum(abs(r0-r1)>EPS)))
     fprintf('r0=r2? diff:%d\n',sum(sum(abs(r0-r1)>EPS)))
     fprintf('r0=r3? diff:%d\n',sum(sum(abs(r0-r3)>EPS)))
     fprintf('r1=r2? diff:%d\n',sum(sum(abs(r1-r2)>EPS)))
     fprintf('r1=r3? diff:%d\n',sum(sum(abs(r1-r3)>EPS)))
     fprintf('r2=r3? diff:%d\n',sum(sum(abs(r2-r3)>EPS)))

     %---compare my vs sys
     Ip=uint8(z3)-uint8(z0);
     h9=figure;imshow(Ip,[]), title([b,'_my_sys_uint8_diff_sig_', num2str(sig)],'Interpreter','none');
     Ip2 = imcomplement(Ip);
     h10=figure;imshow(Ip2,[]), title([b,'_my_sys_uint8_flip_diff_sig_', num2str(sig)],'Interpreter','none');
     Ip3=z3-z0;
     h11=figure;imshow(Ip3,[]), title([b,'_my_sys_double_diff_sig_', num2str(sig)],'Interpreter','none');
     Ip4 = imcomplement(Ip3);
     h12=figure;imshow(Ip4,[]), title([b,'_my_sys_double_flip_diff_sig_', num2str(sig)],'Interpreter','none');

     %----compare zero vs replicate
     Ip5=r0-z0;
     %Ip5=uint8(r0)-uint8(z0);
     h13=figure;imshow(Ip5,[]), title([b,'_diff_sig_', num2str(sig)],'Interpreter','none');
     Ip6 = imcomplement(Ip5);
     h14=figure;imshow(Ip6,[]), title([b,'_diff_flip_sig_', num2str(sig)],'Interpreter','none');

     %----save
     resultsFolder=['results\question3\',b,'\'];
     %imwrite(z0,[resultsFolder, b, '_my_zero_twodConv_sig_', num2str(sig),'.jpg']);
     %imwrite(z1,[resultsFolder, b, '_my_zero_imfilter_sig_', num2str(sig),'.jpg']);
     %imwrite(z2,[resultsFolder,b, '_sys_zero_imfilter_sig_', num2str(sig),'.jpg']);
     %imwrite(z3,[resultsFolder,b, '_sys_zero_imgaussfilt_sig_', num2str(sig),'.jpg']);
     %imwrite(r0,[resultsFolder,b, '_my_replicate_twodConv_sig_', num2str(sig),'.jpg']);
     %imwrite(r1,[resultsFolder,b, '_my_replicate_imfilter_sig_', num2str(sig),'.jpg']);
     %imwrite(r2,[resultsFolder,b, '_sys_replicate_imfilter_sig_', num2str(sig),'.jpg']);
     %imwrite(r3,[resultsFolder, b, '_sys_replicate_imgaussfilt_sig_', num2str(sig),'.jpg']);
     print(h1,'-djpeg',[resultsFolder, b, '_my_zero_twodConv_sig_', num2str(sig),'.jpg']);
     print(h2,'-djpeg',[resultsFolder, b, '_my_zero_imfilter_sig_', num2str(sig),'.jpg']);
     print(h3,'-djpeg',[resultsFolder,b, '_sys_zero_imfilter_sig_', num2str(sig),'.jpg']);
     print(h4,'-djpeg',[resultsFolder,b, '_sys_zero_imgaussfilt_sig_', num2str(sig),'.jpg']);
     print(h5,'-djpeg',[resultsFolder,b, '_my_replicate_twodConv_sig_', num2str(sig),'.jpg']);
     print(h6,'-djpeg',[resultsFolder,b, '_my_replicate_imfilter_sig_', num2str(sig),'.jpg']);
     print(h7,'-djpeg',[resultsFolder,b, '_sys_replicate_imfilter_sig_', num2str(sig),'.jpg']);
     print(h8,'-djpeg',[resultsFolder, b, '_sys_replicate_imgaussfilt_sig_', num2str(sig),'.jpg']);
     print(h9,'-djpeg',[resultsFolder, b, '_my_sys_uint8_diff_sig_', num2str(sig),'.jpg']);
     print(h10,'-djpeg',[resultsFolder, b, '_my_sys_uint8_flip_diff_sig_', num2str(sig),'.jpg']);
     print(h11,'-djpeg',[resultsFolder, b, '_my_sys_double_diff_sig_', num2str(sig),'.jpg']);
     print(h12,'-djpeg',[resultsFolder, b, '_my_sys_double_flip_diff_sig_', num2str(sig),'.jpg']);
     print(h13,'-djpeg',[resultsFolder, b, '_diff_sig_', num2str(sig),'.jpg']);
     print(h14,'-djpeg',[resultsFolder, b, '_diff_flip_sig_', num2str(sig),'.jpg']);

  5. 分析
    值得注意的是，imfilter默认的方法是corr相关，而实现的twodConv是卷积，根据网站上说相关和卷积就差了旋转180度，由于高斯核是对称的，所以在这里效果一样。

    自己和系统方法，在uint8数据类型下，差值图为黑色，说明两个方法得到的结果几乎一致。而在double数据类型下，差值图较为明显，这是因为imshow显示时对颜色重新映射了，所以即使差值很小(<1e-9)，但是重新映射的范围也很小，导致颜色明显。

  6. 参考
    高斯滤波原理：
    https://blog.csdn.net/ytang_/article/details/52749497
    系统高斯核使用：
    https://blog.csdn.net/majinlei121/article/details/50469420
    高斯滤波fspecial代码实现：
    http://www.voidcn.com/article/p-wqwnnifq-rp.html
    https://blog.csdn.net/zhangfuliang123/article/details/76150175

傅里叶变换和逆变换

  1. 知识
    分治法实现fft（重要）；
    https://blog.csdn.net/u011529752/article/details/56496603
    python 图像处理（重要）：
    https://blog.csdn.net/Eastmount/article/details/89474405
    知识：
    https://blog.csdn.net/abcjennifer/article/details/7622228
    傅里叶变换知识:
    https://www.cnblogs.com/mq0036/p/4798315.html

  2. 傅里叶变换
    系统函数：

    1
    2
    F = fft(f) % 一维
    F = fft2(f) % 二维

    自己实现：

    1
    F = dft2D(f);

    dft2D.m函数：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    % method 1
    function F = dft2D(f)
    [rows, cols] = size(f);
    F = zeros(rows, cols);
    for i = 1: rows
      F(i,:) = fft(f(i,:));
    end
    for i = 1: cols
      F(:,i) = fft(F(:,i));
    end
    % method 2
    function F = dft2D_2(f)
    F = fft(fft(f).').';

  3. 傅里叶逆变换
    系统函数；

    1
    2
    f = ifft(F) % 一维
    f = ifft2(F) % 二维

    自己实现:

    1
    2
    f = idft2D(F);
    f = idft2D_2(F);

    idft2D.m和idft2D_2.m函数：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    % method 1
    function f = idft2D(F)
    [rows, cols] = size(F);
    f = zeros(rows, cols);
    for i = 1: rows
       f(i,:) = ifft(F(i,:)); % system function
       %f(i,:) = IFFT(F(i,:));  % my function 
    end
    for i = 1: cols
      f(:,i) = ifft(f(:,i)); % system function
      %f(:,i) = IFFT(f(:,i));  % my function
    end

    % method 2
    function f = idft2D_2(F)
     [M, N] = size(F);
     F_conj = zeros(M, N);
     % --- the same as: F_conj = conj(F)
     for i = 1: M
         for j = 1: N
             if F(i,j)== 0
                 F_conj(i,j) = F(i,j);
             else
                 F_conj(i,j) = abs(F(i,j)).^2./F(i,j);
             end
         end
     end
     f_conj_MN = dft2D(F_conj);
     f_conj = f_conj_MN / M / N;
     % ---- the same as: f = conj(f_conj);
      for i = 1: M
         for j = 1: N
             if F(i,j)== 0
                 f(i,j) = f_conj(i,j);
             else
                 f(i,j) = abs(f_conj(i,j)).^2./f_conj(i,j);
             end
         end
      end

    原始傅里叶逆变换IFFT.m实现：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    function f = IFFT(F)
    sz = size(F);
    if sz(1) > 1 && sz(2) > 1 || numel(sz) >2
      f = -1;
      return;
    end
    M = max(sz(:));
    if M == 1
      f = F;
      return;
    end
    f = zeros(sz(1), sz(2));
    for n = 1:M
      for m = 1:M
         w = exp(i*2*pi*(n-1)*(m-1)/M);
         f(n) = f(n) + F(m) * w;
      end
    end
    f = f./M;

  4. 中心化
    在空间域实现中心化

    1
    2
    3
    4
    5
    6
    7
    % before fft
    img_center = zeros(rows, cols);
    for i = 1: rows
      for j = 1: cols
         img_center(i,j) = img(i,j)*(-1)^(i+j-2);
      end
    end

    在频率域实现中心化

    1
    2
    3
    4
    5
    6
    7
    8
    9
    % after fft
    img_dft2D_center_fftshift = fftshift(img_dft2D); % system function
    img_dft2D_center_circshift = circshift(img_dft2D, [floor(rows/2), floor(cols/2)]); % system function
    % 交换四个现象
    img_dft2D_center_my = zeros(rows, cols);
    img_dft2D_center_my(1:rows/2,1:cols/2) = img_dft2D((rows-rows/2+1):rows, (cols-cols/2+1):cols);
    img_dft2D_center_my(1:rows/2,(cols-cols/2+1):cols) = img_dft2D((rows-rows/2+1):rows, 1:cols/2);
    img_dft2D_center_my((rows-rows/2+1):rows, 1:cols/2) = img_dft2D(1:rows/2,(cols-cols/2+1):cols);
    img_dft2D_center_my((rows-rows/2+1):rows, (cols-cols/2+1):cols) = img_dft2D(1:rows/2,1:cols/2);

  5. 对数变换

    1
    FF = log(abs(F)+1);

  6. 验证
    计算ifft生成图像和源图像差值

    1
    2
    3
    4
    5
    % ---- cal diff between img_dft2D_idft2D and img_origin
     Ip1=uint8(img)-uint8(abs(img_idft2D));
     Ip2 = imcomplement(Ip1);
     Ip3= img-abs(img_idft2D);
     Ip4 = imcomplement(Ip3);

    从多个角度比较

    1
    2
    3
    4
    5
    6
    % --- img->   fft:  dft2D vs fft2
    % --- img->dft2D->  ifft:    ifdt2D vs idft2D_2 vs ifft2D
    % --- img->dft2D-> center:   my vs fftshift vs circshift 
    % --- img->center->dft2D->log 流水线1
    % --- img->dft2D->center->log 流水线2
    % --- diff img_dft2D_idft2D vs img_origin

  7. 参考
    circshift介绍：https://blog.csdn.net/xuexiopencv/article/details/51202794
    fftshift等对比：
    http://www.jeepxie.net/article/84206.html
    fft2: https://blog.csdn.net/double12754/article/details/88393626

matlab常用接口

  1. 常用接口
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
im2double(img) %0-255转0-1
im2uint8
conj(X) ==  abs(X).^2./X 注意X不能为0
% compare
isequal(img_dft2D_center_log, img_center_dft2D_log)
EPS = 1e-9;
fprintf('dft2D=fft2? diff:%d\n',sum(sum(abs(img_dft2D-img_fft2)>EPS)))
% read data
foldername = '\';
filename = 'rose512.tif';
[a,b1,c] = fileparts([foldername,filename]);
img1 = imread([foldername,filename]);
% show
h5=figure;
set(gcf,'position',[400,100,800,600]);
subplot(2,3,4);
imshow(Ip2,[]), title(['(d) ',b,'_uint8_flip_diff'],'Interpreter','none');
% save
resultsFolder=['results\',b,'\'];
print(h1,'-djpeg',[resultsFolder, b, '_diff_fft.jpg']);
  1. 数据分布
1
2
3
4
5
6
7
8
ymin = min(y(:))
ymax = max(y(:))
x = linspace(ymin, ymax, 21);
y2 = reshape(y,1,[]);
yy = hist(y2,x,'grouped')
yy = yy/sum(yy);
figure(1);
bar(x,yy);
赏

谢谢你请我吃糖果

微信

扫一扫，分享到微信

微信分享二维码

tag:

    缺失模块。
    1、请确保node版本大于6.2
    2、在博客根目录（注意不是yilia根目录）执行以下命令：
    npm i hexo-generator-json-content --save

    3、在根目录_config.yml里添加配置： 

      jsonContent:
    meta: false
    pages: false
    posts:
      title: true
      date: true
      path: true
      text: false
      raw: false
      content: false
      slug: false
      updated: false
      comments: false
      link: false
      permalink: false
      excerpt: false
      categories: false
      tags: true

很惭愧<br><br>只做了一点微小的工作<br>谢谢大家