func_detect1.m

用于对天文图像进行去噪和探测

%------------------------------------
%作者：韦海萍
%日期：2007-7-19
%文件名：detect.m
%------------------------------------
function  [err, SourceImage, DestImage,line_out,column_out,direct_out,sudu_out]=func_detect1(str,reportoutput,FileFormat)

%---------------------- for test ------------------------------------------------
%{
clc
clear all
chengxustart=clock
str{1}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U57U_Image\U57U_Image_0071.fit';
str{2}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U57U_Image\U57U_Image_0072.fit';
str{3}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U57U_Image\U57U_Image_0073.fit';
str{4}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U57U_Image\U57U_Image_0074.fit';
str{5}='U57U_Image_0074的参数报告0815.dat';
str{6}='U57U_Image_0074';

reportoutput=1;
FileFormat='ieee-le';
%}

%---------------------- start ------------------------------------------------
Height = 1024;
Width = 1024;
file1 = fopen(str{1},'rb');
file2 = fopen(str{2},'rb');
file3 = fopen(str{3},'rb');
file4 = fopen(str{4},'rb');

if file1 ~= -1 && file2~= -1 && file3~= -1 && file4~= -1
    err=0;    
    
    fseek(file1,2880,'bof');
    fseek(file2,2880,'bof');
    fseek(file3,2880,'bof');
    fseek(file4,2880,'bof');

    %--------------前两帧图像相减，运算；---------------------------------------

    i44 = fread(file1,[Width,Height],'uint16',FileFormat);
    i45 = fread(file2,[Width,Height],'uint16',FileFormat);

    fclose(file1);
    fclose(file2);

    i1 = i44';
    i2 = i45';

    clear i44 i45 file1 file2;
    %-----------对要进行做帧差的图像先进行增强-----------------------------------
   i1 = wiener2(i1,[3,3]);%使用wiener滤波进行图像增强
    i2 = wiener2(i2,[3,3]);%使用wiener滤波进行图像增强
    i12=zeros(Height,Width);
    for x=1:Height
        for y=1:Width
            if abs(i1(x,y)-i2(x,y))<700%500
                i12(x,y)=0;
            else
                i12(x,y)=1;
            end
        end
    end
    clear i1 i2 x y;
    figure(1),imshow(i12,[]);
    [io12,num12]=bwlabel(i12,8); %ig2
    clear i12;
    pixelnum12 = zeros(1,num12);

    for i=1:num12
        pixel=0;
        for x = 1:Height
            for y = 1:Width
                if io12(x,y) == i
                    pixel = pixel +1;
                end
            end
        end
        pixelnum12(i) = pixel;
        if pixelnum12(i)<10 %|| len12(i)<=4 || wid12(i)<=4
            for m=1:Height
                for n=1:Width
                    if io12(m,n)==i
                        io12(m,n)=0;
                    end
                end
            end
        end
    end
    clear pixelnum12 pixel m n x y;
    jieguo=zeros(Height,Width);
    for i=1:Height
        for j=1:Width
            if io12(i,j)~=0
                jieguo(i,j)=1;
            end
        end
    end
    clear io12 i j;
    
    figure(2),imshow(jieguo,[]);

    i5=zeros(Height,Width);
    bz=0;
    for x=1:Height-70
        for y=1:70
            if jieguo(x,y)>0
                for m=10:1:70
                    for n=10:1:70
                        if jieguo(x+m,y+n)==0
                            bz=1+bz;
                        else
                            %                         bz=bz-1;
                            i5(x+m,y+n)=jieguo(x+m,y+n);
                        end
                    end
                end
                
                if (bz==61*61 && i5(x,y)==0)
                    i5(x,y)=0;
                else
                    i5(x,y)=jieguo(x,y);
                end
                bz=0;
            else
                i5(x,y)=0;
                bz=0;
            end
        end
        for y=71:Width-70
            if jieguo(x,y)>0
                for m=10:1:70
                    for n=-70:1:70
                        if abs(n)>10
                            if jieguo(x+m,y+n)==0
                                bz=1+bz;
                            else
                                %                         bz=bz-1;
                                i5(x+m,y+n)=jieguo(x+m,y+n);
                            end
                        end
                    end
                end
                
                if (bz==61*120 && i5(x,y)==0)
                    i5(x,y)=0;
                else
                    i5(x,y)=jieguo(x,y);
                end
                bz=0;
            else
                i5(x,y)=0;
                bz=0;
            end
        end
        for y=Width-69:Width
            if jieguo(x,y)>0
                for m=10:1:70
                    for n=10:1:70
                        if jieguo(x+m,y-n)==0
                            bz=1+bz;
                        else
                            %                         bz=bz-1;
                            i5(x+m,y-n)=jieguo(x+m,y-n);
                        end
                    end
                end
               
                if (bz==61*61 && i5(x,y)==0)
                    i5(x,y)=0;
                else
                    i5(x,y)=jieguo(x,y);
                end
                bz=0;
            else
                i5(x,y)=0;
                bz=0;
            end
        end
    end
    clear jieguo x y m n bz;
  

    %----------------- 后两帧图像相减，运算；------------------------------------
    i44 = fread(file3,[Width,Height],'uint16',FileFormat);
    i45 = fread(file4,[Width,Height],'uint16',FileFormat);

    fclose(file3);
    fclose(file4);

    i1 = i44';
    i2 = i45';
    clear i44 i45 file3 file4;

    %--------------- 首先进行高通滤波，去除图像的背景和星下方的白色亮线；---------------------------------
   
    ig1 = zeros(Height,Width);
    ig2 = zeros(Height,Width);
    for y=1:Width
        for x=10:Height-9
            qsjzh=(i1(x-9,y)+i1(x-8,y)+i1(x-7,y)+i1(x-6,y)+i1(x-5,y)+i1(x-4,y)+i1(x-3,y)+i1(x-2,y))/8;%前面8个行系数的均值；
            hsjzh=(i1(x+2,y)+i1(x+3,y)+i1(x+4,y)+i1(x+5,y)+i1(x+6,y)+i1(x+7,y)+i1(x+8,y)+i1(x+9,y))/8;%后面8个行系数的均值；
            yuzhi=max(qsjzh,hsjzh);
            if i1(x,y)<1.3*yuzhi
                ig1(x,y)=0;
            else
                ig1(x,y)=i1(x,y);
            end
        end
    end

    for y=1:Width
        for x=10:Height-9
            qsjzh=(i2(x-9,y)+i2(x-8,y)+i2(x-7,y)+i2(x-6,y)+i2(x-5,y)+i2(x-4,y)+i2(x-3,y)+i2(x-2,y))/8;
            hsjzh=(i2(x+2,y)+i2(x+3,y)+i2(x+4,y)+i2(x+5,y)+i2(x+6,y)+i2(x+7,y)+i2(x+8,y)+i2(x+9,y))/8;
            yuzhi=max(qsjzh,hsjzh);
            if i2(x,y)<1.3*yuzhi
                ig2(x,y)=0;
            else
                ig2(x,y)=i2(x,y);
            end
        end
    end
    clear qsjzh hsjzh yuzhi x y;
    
    %----------------- 统计包括卫星在内的所有星的个数，并计算其位置和大小；--------------------------

    [io2,num2_t]=bwlabel(ig2,8); %ig2

    len2_t = zeros(1,num2_t);
    wid2_t = zeros(1,num2_t);
    pixelnum2_t = zeros(1,num2_t);
    line2_t = zeros(1,num2_t);
    column2_t = zeros(1,num2_t);
    num2=0;
    for i=1:num2_t
        sumx=0;sumy=0;sumi=0;pixel=0; listx=[];listy=[];
        for x = 1:Height
            for y = 1:Width
                if io2(x,y) == i
                    listx = [listx x];
                    listy = [listy y];
                    pixel = pixel +1;
                    sumx = sumx + x*i2(x,y);
                    sumy = sumy + y*i2(x,y);
                    sumi = sumi + i2(x,y);
                end
            end
        end
        if pixel>1
            num2=num2+1;
            len2_t(num2) = max(listx) - min(listx) + 1;
            wid2_t(num2) = max(listy) - min(listy) + 1;
            pixelnum2_t(num2) = pixel;
            line2_t(num2) = sumx/sumi;
            column2_t(num2) = sumy/sumi;
        end
    end
    clear ig2 io2 listx listy pixel sumx sumy sumi;
   
    pixelnum2 = zeros(1,num2);
    line2 = zeros(1,num2);
    column2 = zeros(1,num2);
    len2 = zeros(1,num2);
    wid2 = zeros(1,num2);
    for i = 1 : num2
        pixelnum2(i) = pixelnum2_t(i);
        line2(i) = line2_t(i);
        column2(i) = column2_t(i);
        len2(i) = len2_t(i);
        wid2(i) = wid2_t(i);
    end
    clear pixelnum2_t line2_t column2_t len2_t wid2_t;

    %-------------- 相邻帧图像相减，灰度差小于700的视为背景；-------------------------------------
    i12=zeros(Height,Width);
    for x=1:Height
        for y=1:Width
            if abs(i1(x,y)-i2(x,y))<700%500
                i12(x,y)=0;
            else
                i12(x,y)=1;
            end
        end
    end 
    clear i1 x y;
    
    %------------- 目标面积小于12时，视为噪声点；---------------------------------------------

    [io12,num12]=bwlabel(i12,8);
    clear i12;
    pixelnum12 = zeros(1,num12);

    for i=1:num12
        pixel=0;
        for x = 1:Height
            for y = 1:Width
                if io12(x,y) == i
                    pixel = pixel +1;
                end
            end
        end
        pixelnum12(i) = pixel;
        if pixelnum12(i)<10 %|| len12(i)<=4 || wid12(i)<=4
            for m=1:Height
                for n=1:Width
                    if io12(m,n)==i
                        io12(m,n)=0;
                    end
                end
            end
        end
    end

    jieguo=zeros(Height,Width);
    for i=1:Height
        for j=1:Width
            if io12(i,j)~=0
                jieguo(i,j)=1;
            end
        end
    end
    
    clear io12 pixelnum12 i j x y;

    %--------------- 搜索星左下和右下区域，[-70,-10] [10,70]；-----------------------------------------

    i4=zeros(Height,Width);
    bz=0;
    for x=1:Height-70
        for y=1:70
            if jieguo(x,y)>0
                for m=10:1:70
                    for n=10:1:70
                        if jieguo(x+m,y+n)==0
                            bz=1+bz;
                        else
                            %  bz=bz-1;
                            i4(x+m,y+n)=jieguo(x+m,y+n);
                        end
                    end
                end                
                if (bz==61*61 && i4(x,y)==0)
                    i4(x,y)=0;
                else
                    i4(x,y)=jieguo(x,y);
                end
                bz=0;
            else
                i4(x,y)=0;
                bz=0;
            end
        end
        for y=71:Width-70
            if jieguo(x,y)>0
                for m=10:1:70
                    for n=-70:1:70
                        if abs(n)>10
                            if jieguo(x+m,y+n)==0
                                bz=1+bz;
                            else
                                % bz=bz-1;
                                i4(x+m,y+n)=jieguo(x+m,y+n);
                            end
                        end
                    end
                end                
                if (bz==61*120 && i4(x,y)==0)
                    i4(x,y)=0;
                else
                    i4(x,y)=jieguo(x,y);
                end
                bz=0;
            else
                i4(x,y)=0;
                bz=0;
            end