% clear;
%Acquire sample data for both target and background. Then Use IRELIEF algorithm to get a weight vector and a threshold which are needed in testing phase. 

%Prepare sample data for the IRELIEF algorithm
% PARA=struct('east_min',1653166,'east_max',1655165,'north_min',7367489,...
%     'north_max',7370488,'n_rows',3000,'n_cols',2000,'r',10);
% east_min=1653166;
% east_max=1655165;
% north_min=7367489;
% north_max=7370488;
% n_rows=3000;
% n_cols=2000;
% % reference order the paper used
% refm1=[3 4 5 3 4 5             %(ref 1)
%       4 5 2 4 5 2
%       5 2 3 5 2 3
%       2 3 4 2 3 4];
% 
% % other references
% refm2=[4 5 3 4 5 3             %(ref 2)
%       5 2 4 5 2 4
%       2 3 5 2 3 5
%       3 4 2 3 4 2];
% refm3=[5 3 4 5 3 4               %(ref 3)
%       2 4 5 2 4 5
%       3 5 2 3 5 2
%       4 2 3 4 2 3];
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% r=7;    %Window size = 2*r+1
% thresh_1=0.25;  %pixels with intensity smaller than thresh_1 
%                 %will be set to 0 directly
% %refm=refm1;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dim=(2*r+1)^2;
% nb_targ=2;
% PA=[1653166 1655165 7367489 7370488 3000 2000 r];
% 
% mis=3; pas=5;
% nh=5;
% h=ones(nh,nh)/(nh^2);
% nborder=round((nh-1)/2);  


      
% info = VHF_get_image_info();
% startx=0;
% for mission=2:3
%     for pass=1:6
%         startx=((mission-2)*6+pass-1)*25;
%         fn=sprintf('D:\\LabWork\\PublicData\\VHF Images\\Mission%d\\pass%d',mission,pass);
%         im1 = VHF_read_image(fn,info.n_cols,info.n_rows,1,info.n_cols,1,info.n_rows);
%         fn=sprintf('D:\\LabWork\\PublicData\\VHF Images\\Mission%d\\pass%d',refm(mission-1,pass),pass);
%         ref = VHF_read_image(fn,info.n_cols,info.n_rows,1,info.n_cols,1,info.n_rows);
%         dif=im1-ref;
%         dif=dif.*(dif>0);
%         fn=sprintf('D:\\LabWork\\PublicData\\VHF Images\\Mission%d\\Targets.txt',mission);
%         targets=VHF_read_target_list(fn);
%         minr=n_rows; minc=n_cols;
%         maxr=1; maxc=1;
%         for i=1:targets.N_targets
%             row=round(north_max-targets.N_coord(i)+1);
%             col=round(targets.E_coord(i)-east_min+1);
%             if row<minr
%                 minr=row;
%             end
%             if row>maxr
%                 maxr=row;
%             end
%             if col<minc
%                 minc=col;
%             end
%             if col>maxc
%                 maxc=col;
%             end
%             target{i}=dif(row-r:row+r,col-r:col+r);
%             x(startx+i,:)=reshape(target{i},1,dim);
%         %     figure,imshow(target{i});
%         end
%     end
% end
% x1=x';  %target samples


%Use difference image generated from mission 3 pass 5 with mission x pass 5
%as reference. Target positions are from Karl.

% fid=fopen('Difference\m3p5_m5p5.bin','r');
% di=fread(fid,[n_rows n_cols], 'single');
% dif=di.*(di>0);
% nh=5;
% h=ones(nh,nh)/(nh^2);
% dif_conv=conv2(dif,h);
% nborder=round((nh-1)/2);
% %denoise
% %dif=single(dif_conv(nborder+1:size(dif_conv,1)-nborder,nborder+1:size(dif_conv,2)-nborder));
% dif=single(dif);
% fclose(fid);
% clear di




fn=sprintf('Mission%d\\pass%d',mis,pas);
info = VHF_get_image_info();
im1 = VHF_read_image(fn,info.n_cols,info.n_rows,1,info.n_cols,1,info.n_rows);
fn=sprintf('Mission%d\\pass%d',refm(mis-1,pas),pas);
im2 = VHF_read_image(fn,info.n_cols,info.n_rows,1,info.n_cols,1,info.n_rows);
fn=sprintf('Mission%d\\Targets.txt',mis);
deploy=VHF_read_target_list(fn);

%Get the difference
dif=im1-im2;
dif=dif.*(dif>thresh_1);  %only consider the positive part
% dif=single(dif);
dif=conv2(dif,h);
dif=(dif(nborder+1:size(dif,1)-nborder,nborder+1:size(dif,2)-nborder));                
clear im1 im2




% fn='target_lists\Adolf_Fredrik.Targets.txt';
% Adolf=VHF_read_target_list(fn);
% fn='target_lists\Fredrik.Targets.txt';
% Fredrik=VHF_read_target_list(fn);
fn='target_lists\Karl.Targets.txt';
Karl=VHF_read_target_list(fn);
% fn='target_lists\Sigismund.Targets.txt';
% Sigismund=VHF_read_target_list(fn);

minr=n_rows; minc=n_cols;
maxr=1; maxc=1;
stx=1;
clear x1;
for i=1:deploy.N_targets
    row=round(north_max-deploy.N_coord(i)+1);
    col=round(deploy.E_coord(i)-east_min+1);
%     if row<minr
%         minr=row;
%     end
%     if row>maxr
%         maxr=row;
%     end
%     if col<minc
%         minc=col;
%     end
%     if col>maxc
%         maxc=col;
%     end
%     target{i}=dif(row-r:row+r,col-r:col+r);
    for ii=row-nb_targ:row+nb_targ
        for jj=col-nb_targ:col+nb_targ
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Rotation invariant tranform
            p=dif(ii-r:ii+r,jj-r:jj+r);
            Rp=radon(p,theta);
            FRp=abs(fft(Rp,[],2));
            p=reshape(FRp,1,numel(FRp));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
%             p=reshape(dif(ii-r:ii+r,jj-r:jj+r),1,dim);
%            Do_curvelet;
%             Do_FFT;
            x1(stx,:)=p;
%             p=[p1 p2];
%             if [1 p]*surface>1.5
%                 x(stx,1:size(p1,2))=p1;
%                 x(stx,size(p1,2)+1:size(p1,2)+size(p2,2))=p2;
                stx=stx+1;
%             end
        end
    end
%     figure,imshow(target{i});
end
x1=x1';  %target samples

%Background positions are manually marked as following
% bgx=[223 480 804 1010 1644 1372 240 821 1769 809 268 268 690 869 665 291 139 1499 994 673 576 1154 812 864 1403];
% bgy=[216 208 308 239 308 820 1112 1156 2716 2482 2724 2393 2056 1603 1204 892 690 1476 1813 1685 1348 2201 2496 1777 2005];
% 
% nbg=size(bgx,2);
% for i=1:nbg
%     x2(:,i)=reshape(dif(bgy(i)-r:bgy(i)+r,bgx(i)-r:bgx(i)+r),dim,1); %background samples
% end

%Get background samples, all pixels besides targets.
border=100;
idx=1;
clear x2;
for ii=650:40:1650%r+1:20:n_rows-r
    for jj=r+1:40:n_cols-r
%         if i>minr-border && i<maxr+border && j>minc-border && j<maxc+border
%             continue; % skip the real target region
%         end
%         p=reshape(dif(ii-r:ii+r,jj-r:jj+r),1,dim);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Rotation invariant tranform
            p=dif(ii-r:ii+r,jj-r:jj+r);
            Rp=radon(p,theta);
            FRp=abs(fft(Rp,[],2));
            p=reshape(FRp,1,numel(FRp));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%        Do_curvelet;
%         Do_FFT;
        x2(idx,:)=p;
%         x2(idx,:)=p;
%         if [1 p]*surface>1.5
%             x2(idx,:)=p;
%         end
        idx=idx+1;
    end
end
x2=x2';

xx=[x1 x2]; %sample data for IRELIEF algorithm
targ=[ones(1,size(x1,2)) 2*ones(1,size(x2,2))]; %labels of the sample data

%Set parameters for IRELIEF algorithm
para.it=300;
para.distance='euclidean';
para.kernel='exp';
para.sigma=25;
para.Outlier=1;
[WEIGHT, history] = IRELIEF(xx, targ, para);
%y1, y2 are projected 1-D samples
y1=WEIGHT'*x1;
y2=WEIGHT'*x2;
plot(y1,zeros(1,size(y1,2)),'bx',y2,ones(1,size(y2,2)),'ro');
xtarg=mean(x1')'.*WEIGHT;
xback=mean(x2')'.*WEIGHT;
