大气湍流退化图像复原技术研究及DSP实现

短命女 2024-04-01 09:53 54阅读 0赞

第一章 绪论 1  
1.1 研究背景 1  
1.2 国内外研究现状 1  
1.3 本文的研究内容 5  
1.4 本文的组织结构 5  
第二章 大气湍流退化图像复原技术 7  
2.1 图像退化及复原数学模型概述 7  
2.1.1 图像退化的数学模型 7  
2.1.2 图像复原的数学模型 8  
2.2 大气湍流退化图像概述 9  
2.2.1 大气湍流的基本物理特性 9  
2.2.2 长曝光大气湍流退化图像 10  
2.2.3 短曝光大气湍流退化图像 11  
2.3 主流大气湍流退化图像复原算法 13  
2.3.1 RL-IBD 算法 13  
2.3.2Lucky Region 算法 15  
2.3.3Speckle Image 算法 18  
2.3.4 APEX 算法 22  
2.4 小结 25  
第三章 基于图像频谱特征的大气湍流退化点扩散函数估计 26  
3.1 点扩散函数估计的一般方法 26  
3.2 基于自然景物频谱特点的大气湍流点扩散函数估计方法 27  
3.2.1 自然景物图像频谱分析 27  
3.2.2 近似等腰三角形频谱重建模型建立 38  
3.2.3 大气湍流退化点扩散函数估计流程 42  
3.3 小结 43  
第四章 大气湍流退化图像复原及实验结果分析 44  
4.1 大气湍流退化图像复原 44  
4.1.1 大气湍流退化图像复原过程归纳 44  
4.1.2 大气湍流退化图像复原实验 44  
4.2 大气湍流退化图像复原结果客观评价 49  
4.2.1 图像复原客观评价指标 49  
4.2.2 大气湍流退化图像复原客观评价计算结果 50  
4.3 小结 50  
第五章 大气湍流退化图像复原技术 DSP 实现 52  
5.1 硬件需求分析 52  
5.2 硬件构架设设计 53  
5.3 PCB 设计 55  
5.3.1 电路原理图设计 55  
5.3.2 PCB 设计 56  
5.4 程序设计以及软硬件调试 56  
5.4.1 硬件调试 56  
5.4.2 软件设计及调试 57  
5.5 硬件实现结果 59  
5.6 小结 59  
第六章 总结与展望 61  
参考文献 63  
致谢 69

1.3 本文的研究内容

依据本课题的应用背景，本文的研究对象为大视场、远距离条件下拍摄的因大气湍 流而退化的自然图像。本文首先研究了图像退化和复原的数学模型，然后调研了当前主 流的大气湍流退化图像复原算法，分析了各自的优缺点。其次研究了大气湍流导致光学 系统成像退化的物理过程和基本的数学关系，为大气湍流退化点扩散函数的建立奠定了 一定的基础。随后针对大视场、远距离条件下拍摄的自然图像，本文通过数学模拟其退 化过程，对比分析其频谱的变化规律，并得出了能用于点扩散函数估计的重要特征。利 用这些特征本文提出了大视场、远距离条件下拍摄的清晰自然图像频谱的近似等腰三角 形简化模型，并利用该模型结合退化图像的频谱信息估计出了大气湍流退化图像的点扩 散函数，最后利用该点扩散函复原出了清晰图像。本文还在基于DSP+FPGA的数字图像 实时处理硬件系统方面进行了探究，并开发出了一套硬件系统，同时在该硬件系统中验 证了本文提出的图像复原算法。  
1.4 本文的组织结构  
本文按如下的结构进行展开： 第一部分为绪论部分，主要介绍了本课题研究背景和意义以及国内外研究现状，并  
交代了本文的组织结构。 第二部分完成研究的前导工作，分析需求与研究的难点，找到比较切合需求的已有  
算法进行实现与测试，分析已有算法在应用中的优劣势。 第三部分在前导工作的基础上，分析了大量的自然图像的频谱特征，找出传统点扩  
散函数估计算法在本应用中的难点和关键问题，并提出符合本应用需求的新的点扩散函 数估计算法，给出了设计思路与实现步骤。  
第四部分介绍了大气湍流图像复原在MATLAB软件中的实验情况，详细说明了实验 条件，给出了实验结果，设计了客观评价标准，并依据该标准分析了图像复原前后质量 提高的程度，以证实本文算法的有效性。  
第五部分介绍了本课题开发的硬件处理系统，详细介绍了系统电路原理图PCB以及 软硬件调试的内容和关键点，最后给出了本文算法在硬件中的实现的结果。  
第六部为全文总结，并针对现有研究仍存在的问题给出了改进建议，对后续工作情况进行了一定的规划。  
5.1 硬件需求分析

本硬件系统需要根据上位机发送的指令要求完成对标准视频信号的采集（将模拟视频 信号进行数字化）、存储、处理（完成图像复原）以及输出（将数字图像数据转化为模 拟视频信号），且整个过程必须满足实时性的要求。  
系统数据量分析：  
该系统的输入信号为标准 PAL 制视频信号，其信号格式符合国际电信协会提出的  
ITU-R601/656 （CCIR601/656）中的 BT.656 标准。BT.656 除了传输 4:2:2 的 YCbCr 视  
频数据流外，本文转载自http://www.biyezuopin.vip/onews.asp?id=12691还有行、列同步所用的控制信号\[51\]。如图 34 所示，一帧图像数据由一个  
625 行、每行 1728 字节的数据块组成。其中，23～311 行是偶数场视频数据，336～624 行是奇数场视频数据，其余为垂直控制信号。  
图 35 中，每行数据包含水平控制信号和 YCbCr 视频数据信号。视频数据信号排列顺序 为 Cb-Y-Cr-Y。每行开始的 288 字节为行控制信号，开始的 4 字节为 EAV 信号(有效视频结束)，紧接着 280 个固定填充数据，最后是 4 字节的 SAV 信号(有效视频起始)。SAV 和 EAV 信号有 3 字节的前导：FF、FF、00；最后 1 字节 XY 表示该行位于整个数据帧 的位置及如何区分 SAV、EAV。XY 字节各比特位含义如图 36。

clear all;close all;
    Image_Input=imread('19.bmp');
    dims=ndims(Image_Input);                     %获取矩阵的维数
    if dims==3
       Image_Gray=rgb2gray(Image_Input);           
    elseif dims~=2
        return;
    end
    [h,w]=size(Image_Gray);                  
    xr=8;
    xl=8;
    yt=8;
    yb=8;
    width=w-xr-xl;
    height=h-yt-yb;
    Image_Cut=imcrop(Image_Gray,[xr,yt,width-1,height-1]);   %  切除黑边
    figure(1);
    imshow(Image_Cut,[]);
    [Row,Column]=size(Image_Cut);
    Image_fft=fftshift(fft2(Image_Cut));    %傅里叶变换,平移
    Image_abs=abs(Image_fft);    %求频谱
    Image_norm=Image_abs/max(max(Image_abs));     %频谱归一化
    Image_ln=log10(Image_norm);     %对数缩放处理
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Image_array=Image_ln(fix(Row/2),1:end);      %获取水平中心处平频谱
    % if rem(Column,2)==0                            %使用中心零点对称坐标显示
    %     X=-fix(Column/2)+1:1:fix(Column/2);
    % else
    %     X=-fix(Column/2):1:fix(Column/2);
    % end
    % figure(3);
    % plot(X,Image_array);
    % xlabel('离散频率')
    % ylabel('傅里叶频率归一化后取自然对数');
    % Datablur=Image_array(fix(Column/2):end);
    % figure(4);
    % semilogx(Datablur,'b:*');
    % xlabel('离散频率的自然对数');
    % ylabel('傅里叶频率归一化后取自然对数');
    % hold on;
    % X_half=fix(Column/2):Column;
    % count=size(X_half);
    % k=1:count(2);
    % p=polyfit(log10(k),Datablur,1);
    % line=polyval(p,log10(k));
    % semilogx(line,'r:*');
    % hold off;
    if rem(Row,2)~=0         %获取水平中心处平频谱
        Image_array=Image_ln((Row+1)/2,1:end);      
    else
        Image_array=(Image_ln(Row/2,1:end)+Image_ln(Row/2+1,1:end))/2;
    end
    if rem(Column,2)==0                            %使用中心零点对称坐标显示
        Data_X=-Column/2:Column/2-1;
    else
        Data_X=-(Column-1)/2:(Column-1)/2;
    end
    figure(2);
    plot(Data_X,Image_array);
    xlabel('离散频率')
    ylabel('傅里叶频率归一化后取自然对数');
    if rem(Column,2)==0
        Data_Y=Image_array(Column/2+1:end);
    else
        Data_Y=Image_array((Column+1)/2:end);
    end
    figure(3);
    semilogx(Data_Y,'b:*');
    xlabel('离散频率的自然对数');
    ylabel('傅里叶频率归一化后取自然对数');
    hold on;
    if rem(Column,2)==0
        k=1:Column/2;
    else
        k=1:(Column+1)/2;
    end
    p=polyfit(log10(k),Data_Y,1);
    line=polyval(p,log10(k));
    semilogx(line,'r:*');
    hold off;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 直接将双对数坐标下的拟合直线斜率设为1  重建清晰图像的频谱  并计算参数
    slope=-1;      % 给定双对数最佳拟合曲线的斜率，为大量统计得到的经验值
    focal_point=7;  % 表示清晰图像频谱与退化图像频谱在双对数坐标系的下焦点力原点的距离
    % 假设清晰图像频谱在双对数坐标系下的直线方程为 Y=aX+b,a是斜率，b是截距（待求解）
    % Y=ln(G(0,focal_point)),X=ln(foacl_point),b=Y-aX
    a=slope;
    if rem(Column,2)==0
        x_length=Column/2;
        Y=Image_array(Column/2+focal_point);
        NewSpectrumHalf=zeros(1,Column/2);
    else x_length=(Column-1)/2;
         Y=Image_array((Column+1)/2+focal_point);
         NewSpectrumHalf=zeros(1,(Column-1)/2);
    end
    X=log10(focal_point);
    b=Y-slope*X;
    for i=1:x_length
        x_temp=log10(i);
        if i<focal_point
            if rem(Column,2)==0
                NewSpectrumHalf(i)=Image_array(Column/2+i); 
            else NewSpectrumHalf(i)=Image_array((Column+1)/2+i);
            end
        else
            NewSpectrumHalf(i)=a*x_temp+b;
        end
    end
    New_spectrum_1=zeros(1,Column);
    for i=1:Column                           
        if rem(Column,2)==0
            if i<=Column/2
                New_spectrum_1(i)=NewSpectrumHalf(1+Column/2-i);
            else
                New_spectrum_1(i)=NewSpectrumHalf(i-Column/2);
            end
        else
            if i<(Column+1)/2
                New_spectrum_1(i)=NewSpectrumHalf((1+Column)/2-i);
            elseif i==(Column+1)/2
                New_spectrum_1(i)=Image_array(i);
            else
                New_spectrum_1(i)=NewSpectrumHalf(i-(Column+1)/2);
            end
        end
    end
    figure(4);   
    plot(Data_X,Image_array);
    hold on;
    plot(Data_X,New_spectrum_1,'r:*');
    hold off;
    Dvalue_1=Image_array-New_spectrum_1;
    figure(5);
    plot(Data_X,Dvalue_1);
    xlabel('离散频率'); 
    ylabel('含参数的代数式');
    scale=100;
    Data_X_stract=[fix(Column/2)-scale:fix(Column/2)+scale]-fix(Column/2);
    Dvalue_stract_1=Dvalue_1(Data_X_stract+fix(Column/2));
    figure(6);
    plot(Data_X_stract,Dvalue_stract_1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 用直线重建清晰图像的频谱y=ax+b  并计算参数
    Number1=10;                                                  
    Number2=20;                                               
    Number3=5;                                                   % Number3用于修正直线的斜率
    y1temp=Image_array(fix(Column/2)-Number1:fix(Column/2));     
    y_1=sum(y1temp)/(Number1+1);                                 % 注意这里的数据个数不要搞错了，是Number1+1个
    x1=fix(Column/2)-Number3;
    x11=fix(Column/2)+1+Number3;
    y2temp=Image_array(1:Number2);
    y_2=sum(y2temp)/Number2+0.6;                                 % 若采用单张图像进行估计需要将末端整体上台一个1   因为退化前后末端的频谱不可能还是一样。                            
    x2=1;                                                        
    x22=Column;
    coeff_a=(y_1-y_2)/(x1-x2);                     % 对称轴左边
    coeff_b=((y_1+y_2)-coeff_a*(x1+x2))/2;
    coeff_c=(y_1-y_2)/(x11-x22);                   % 对称轴右边
    coeff_d=((y_1+y_2)-coeff_a*(x11+x22))/2;
    for i=1:Column
        if i<=x1
             New_spectrum_2(i)=coeff_a*i+coeff_b;
        elseif i>x1 && i<x11
             New_spectrum_2(i)=Image_array(i);
        elseif i>=x11
             New_spectrum_2(i)=coeff_a*(Column-i)+coeff_b;
        end
    end
    figure(7);
    plot(Data_X,Image_array);
    hold on;
    plot(Data_X,New_spectrum_2,'r*');
    xlabel('离散频率');                             
    ylabel('图像频谱归一化的自然对数');
    hold off
    Dvalue_2=Image_array-New_spectrum_2;
    figure(8);
    plot(Data_X,Dvalue_2);
    xlabel('离散频率'); 
    ylabel('含参数的代数式');
    Dvalue_stract_2=Dvalue_2(Data_X_stract+fix(Column/2));
    figure(9);
    plot(Data_X_stract,Dvalue_stract_2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % alpha=0.001;
    alpha=0.0008195 ;
    belta=5/6;
    Image_deblur=ones(Row,Column);
    if rem(Row,2)==0
        for u=-Row/2:Row/2-1
            if rem(Column,2)==0
                for v=-Column/2:Column/2-1
                    H_Estimate(u+Row/2+1,v+Column/2+1)=exp(-alpha*(u^2+v^2)^belta);
                end
            else 
                for v=-(Column-1)/2:(Column-1)/2
                    H_Estimate(u+Row/2+1,v+(Column-1)/2+1)=exp(-alpha*(u^2+v^2)^belta);
                end
            end
        end
    else
        for u=-(Row-1)/2:(Row-1)/2
            if rem(Column,2)==0
                for v=-Column/2:Column/2-1
                    H_Estimate(u+(Row-1)/2+1,v+Column/2+1)=exp(-alpha*(u^2+v^2)^belta);
                end
            else
                for v=-(Column-1)/2:(Column-1)/2
                    H_Estimate(u+(Row-1)/2+1,v+(Column-1)/2+1)=exp(-alpha*(u^2+v^2)^belta);
                end
            end
        end
    end
    % Image_deblur=Degeneration./H_Estimate;
    % 
    % Image_deblur=ifft2(ifftshift(Image_deblur));
    % Image_deblur=real(Image_deblur);
    % Image_deblur=0.5*(abs(Image_deblur)+Image_deblur);
    % figure(11);
    % imshow(Image_deblur,[]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    h_estimate=real(ifft2(ifftshift(H_Estimate))); % 转化到空域上来
    result=deconvwnr(Image_Cut,h_estimate,0.001);
    result=ifftshift(result); % 对图像进行1、3象限对调，2、4象限对调
    figure(10);
    imshow(result,[]);
    title('维纳滤波复原图像');

![在这里插入图片描述][a06074a9c7b940f48320c11e55f064d4.png_pic_center]  
![在这里插入图片描述][f94075df5d5e4c53b4c2ba04174c05b8.png_pic_center]  
![在这里插入图片描述][397e1a2fd619494a999f3b611d3960cd.png_pic_center]  
![在这里插入图片描述][676059c56d6c41c39d74dc00c495ec82.png_pic_center]  
![在这里插入图片描述][a77a0617d177474eb371b653f45c135b.png_pic_center]  
![在这里插入图片描述][ae58a50c6bcd4ef5b84d8fe3c52843f7.png_pic_center]  
![在这里插入图片描述][5d2ac0b2dba84f96b55ec6cb9f3b4c35.png_pic_center]  
![在这里插入图片描述][fb0b1963d6204e9cad4260f24dd435be.png_pic_center]  
![在这里插入图片描述][a1686aa595b847d5ab108bb17934a04a.png_pic_center]  
![在这里插入图片描述][2c9c1163ca33492992a77a2de22ad152.png_pic_center]

[a06074a9c7b940f48320c11e55f064d4.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/8bf5bdb4aa7f44718f6c83bab9c98936.png
[f94075df5d5e4c53b4c2ba04174c05b8.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/f3948530eed442c8a4777a27b0d12b9f.png
[397e1a2fd619494a999f3b611d3960cd.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/398af6b4c84c4a9ba3b0ae4a1cde9201.png
[676059c56d6c41c39d74dc00c495ec82.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/d3626d1b7e85407596a11b26e35ac45e.png
[a77a0617d177474eb371b653f45c135b.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/f59f242041e2490b97432b31839d69af.png
[ae58a50c6bcd4ef5b84d8fe3c52843f7.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/b2a73b12e44a44b2a589fa83e7b1fe00.png
[5d2ac0b2dba84f96b55ec6cb9f3b4c35.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/40ab5833359347b6a1082dfc387e8cf7.png
[fb0b1963d6204e9cad4260f24dd435be.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/d94a9927001e4c3b88ef190475c67cf8.png
[a1686aa595b847d5ab108bb17934a04a.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/f2b0bae0826f414488b1e223e246d8da.png
[2c9c1163ca33492992a77a2de22ad152.png_pic_center]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/01/a50aeacdf0ab49b3994a9b823640e9c1.png