func_detect2.m

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

%------------------------------------
%作者：韦海萍
%日期：2007-7-19
%文件名：detect.m
% %------------------------------------
function  [err, SourceImage,DestImage,line_out,column_out,direct_out,sudu_out]=func_detect2(str,reportoutput,line_in,column_in,direct_in,sudu_in,FileFormat)
% clc
% clear all
% str{3}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U56U_Image\U56U_Image_0052.fit';
% str{4}='G:\projects\matlab works\863-紫金山天文台星图处理070531\U56U_Image\U56U_Image_0053.fit';
% str{5}='U56U_Image_0053的参数报告0815.dat';
% str{6}='U56U_Image_0053';
% 
% line_in=527.8;
% column_in=182;
% direct_in=111.43;
% sudu_in=35.053;
% 
% reportoutput=2;

Height = 1024;
Width = 1024;
file1 = fopen(str{3},'rb');
file2 = fopen(str{4},'rb');

if file1 ~= -1 && file2~= -1
    err=0;       
    fseek(file1,2880,'bof');
    fseek(file2,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;
    
    ig2 = zeros(Height,Width);
    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);  %ig2(x,y)
                    sumy = sumy + y*i2(x,y);  %ig2(x,y)
                    sumi = sumi + i2(x,y);  %ig2(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的视为背景；-------------------------------------
    jieguo=zeros(Height,Width);
    for x=1:Width
        for y=1:Height
            if abs(i1(x,y)-i2(x,y))<700%500
                jieguo(x,y)=0;
            else
                jieguo(x,y)=1;
            end
        end
    end    

    MoBan=zeros(Height,Width);
    CanKaoNum=size(line_in,2);
    valid=[];    
    weizhi=[];
    geshu=0;
    deltax=zeros(1,CanKaoNum);
    deltay=zeros(1,CanKaoNum);
    NewLine=zeros(1,CanKaoNum);
    NewColumn=zeros(1,CanKaoNum);
    for i=1:CanKaoNum
        deltax(i)=-sudu_in(i)*sind(direct_in(i));
        deltay(i)=sudu_in(i)*cosd(direct_in(i));
        NewLine(i)=line_in(i)+deltax(i);
        NewColumn(i)=column_in(i)+deltay(i);
        valid_cnt=0;
        if NewLine(i)>0 && NewLine(i)<=Height && NewColumn(i)>0 && NewColumn(i)<=Width     
            NewLine(i)=round(NewLine(i));
            NewColumn(i)=round(NewColumn(i));  
            if NewLine(i)>=4 && NewLine(i)<=Height-3 && NewColumn(i)>=4 &&  NewColumn(i)<=Width-3
                for xx = -3:1:3
                    for yy = -3:1:3
                        if jieguo(NewLine(i)+xx,NewColumn(i)+yy) == 1
                            MoBan(NewLine(i)+xx,NewColumn(i)+yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)<4 && NewColumn(i)>=4 && NewColumn(i)<=Width-3
                for xx = 1:1:7
                    for yy = -3:1:3
                        if jieguo(xx,NewColumn(i)+yy) == 1
                            MoBan(xx,NewColumn(i)+yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)>=4 && NewLine(i)<=Height-3 && NewColumn(i)<4
                for xx = -3:1:3
                    for yy = 1:1:7
                        if jieguo(NewLine(i)+xx,yy) == 1
                            MoBan(NewLine(i)+xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)>Height-3 && NewColumn(i)>=4 &&  NewColumn(i)<=Width-3
                for xx = Height-6:1:Height
                    for yy = -3:1:3
                        if jieguo(xx,NewColumn(i)+yy) == 1
                            MoBan(xx,NewColumn(i)+yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)>=4 && NewLine(i)<=Height-3 && NewColumn(i)>Width-3
                for xx = -3:1:3
                    for yy = Width-6:1:Width
                        if jieguo(NewLine(i)+xx,yy) == 1
                            MoBan(NewLine(i)+xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)<4 && NewColumn(i)<4
                for xx = 1:1:7
                    for yy = 1:1:7
                        if jieguo(xx,yy) == 1
                            MoBan(xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)<4 && NewColumn(i)>Width-3
                for xx = 1:1:7
                    for yy = Width-6:1:Width
                        if jieguo(xx,yy) == 1
                            MoBan(xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            elseif NewLine(i)>Height-3 && NewColumn(i)<4
                for xx = Height-6:1:Height
                    for yy = 1:1:7
                        if jieguo(xx,yy) == 1
                            MoBan(xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            else
                for xx = Height-6:1:Height
                    for yy = Width-6:1:Width
                        if jieguo(xx,yy) == 1
                            MoBan(xx,yy) = i;
                            valid_cnt=valid_cnt+1;
                        end
                    end
                end
            end            
            if valid_cnt ~=0
                weizhi = [weizhi i];
                valid = [valid valid_cnt];
                geshu=geshu+1;
            end                
        end
    end
    clear deltax deltay NewLine NewColumn;
%     for i=1:Height
%         for j=1:Width
%             temp(i,j)=MoBan(i,j) & jieguo(i,j);
%         end
%     end
    direct = zeros(1,num2);
    sudu = zeros(1,num2);
    
    if geshu==0
        line_out=[];
        column_out=[];
        direct_out=[];
        sudu_out=[];
    else
        line_out=zeros(1,geshu);
        column_out=zeros(1,geshu);
        deltax=zeros(1,geshu);
        deltay=zeros(1,geshu);
        jiaodu=zeros(1,geshu);
        direct_out=zeros(1,geshu);
        sudu_out=zeros(1,geshu);
        for i=1:geshu            
            sumx=0;sumy=0;sumi=0;
            for x = 1:Height
                for y = 1:Width
                    if MoBan(x,y) == weizhi(i)
                        sumx = sumx + x*i2(x,y);
                        sumy = sumy + y*i2(x,y);
                        sumi = sumi + i2(x,y);
                    end
                end
            end
            line_out(i) = sumx/sumi;
            column_out(i) = sumy/sumi;

            deltax(i) = line_out(i)-line_in(weizhi(i));
            deltay(i) = column_out(i)-column_in(weizhi(i));
            jiaodu(i) = atand(abs(deltax(i))/abs(deltay(i)));

            if deltax(i)>0 && deltay(i)>0
                direct_out(i) = 360-jiaodu(i);
            elseif deltax(i)>0 && deltay(i)<0
                direct_out(i) = 180+jiaodu(i);
            elseif deltax(i)<0 && deltay(i)>0
                direct_out(i) = jiaodu(i);
            else
                direct_out(i) = 180-jiaodu(i);
            end
            sudu_out(i) = sqrt(deltax(i).^2+deltay(i).^2);
        end
        for i=1:num2
            for j=1:geshu
                if abs(round(line2(i))-round(line_out(j))) <= 2 && abs(round(column2(i))-round(column_out(j))) <= 2
                    direct(i)=direct_out(j);
                    sudu(i)=sudu_out(j);
                end
            end
        end
    end  
    
    SourceImage=uint16(i2);
    DestImage=zeros(Height,Width);
    for i=1:Height
        for j=1:Width
            if MoBan(i,j) ~=0
                DestImage(i,j)=1;
            end
        end
    end

    if reportoutput == 2
        %--------- 写参数报告dat文件 ---------------------------------------
        fid1=fopen(str{5},'w');
        filename=strcat(str{6},'.fit');
        fprintf(fid1,'FILENAME=%s \n',filename);
        fseek(fid1,0,'cof');
        fprintf(fid1,'TARGET=%s \n','Image');
        fseek(fid1,0,'cof');
        fprintf(fid1,'HEIGHT=%d \n',Height);
        fseek(fid1,0,'cof');
        fprintf(fid1,'WIDTH=%d \n',Width);
        fseek(fid1,0,'cof');
        fprintf(fid1,'STAROBJ_NUM=%d \n',num2);
        fseek(fid1,0,'cof');
        fprintf(fid1,'SATEOBJ_NUM=%d \n',geshu);
        fseek(fid1,0,'cof');
        mubiao=[];
        for i = 1:num2
            if sudu(i) == 0
                fprintf(fid1,'Line=%10.4f   Column=%10.4f   Pixelnum=%6d  Star_Len=%6d  Star_Wid=%6d  Direct=%6.2f   V=%6.2f \n',line2(i),column2(i),pixelnum2(i),len2(i),wid2(i),direct(i),sudu(i)); %
                fseek(fid1,0,'cof');
            else
                mubiao = [mubiao i];
            end
        end
        %-------- 将卫星的信息放在后面输出 -------------------------
        if geshu>0
            %     if geshu==1
            %         fprintf(1,'There is %d target.\n',geshu);
            %     else
            %         fprintf(1,'There are %d targets.\n',geshu);
            %     end
            for i=1:geshu
                fprintf(fid1,'Line=%10.4f   Column=%10.4f   Pixelnum=%6d  Star_Len=%6d  Star_Wid=%6d  Direct=%6.2f   V=%6.2f \n',line2(mubiao(i)),column2(mubiao(i)),pixelnum2(mubiao(i)),len2(mubiao(i)),wid2(mubiao(i)),direct(mubiao(i)),sudu(mubiao(i))); %
                fseek(fid1,0,'cof');
                %         fprintf(1,'The information of target %d are:\n',i);
                %         fprintf(1, ' target_line=%f ;\n target_column=%f ;\n target_size=%d ;\n target_len=%d ;\n target_wid=%d ;\n target_direct=%f ;\n target_speed=%f.\n',line2(mubiao(i)),column2(mubiao(i)),pixelnum2(mubiao(i)),len2(mubiao(i)),wid2(mubiao(i)),direct(mubiao(i)),sudu(mubiao(i)));
            end
            % else
            %     fprintf(1,'There is no target.\n');
        end
        fclose(fid1);
    end
else
    err=1;
    SourceImage=zeros(Height,Width);
    DestImage=zeros(Height,Width);
    line_out=[];
    column_out=[];
    direct_out=[];
    sudu_out=[];
end