clear;
close all;

%{
Suppose a matrix "fir" represents the population responses, with each row representing a stimulus, and each column representing a cell.
Stimuli 1-800: similarly-sized fragment; 801-810: gray screen; 811-910: varying-sized fragment; 911-1140: complete face; 1141-1940: occluded face
Leftright indicates in which hemisphere the neuron was recorded in: 0-right; 1-left
The 6-D facial features of each stimulus were saved as a matrix "aps", where each row represents a feature dimension and each column represents a stimulus.
The location and size of the fragment or the occluder were saved as a matrix "asu"; three rows indicate the x,y locations, and the radii, respectively; each row represents one stimulus.
%}

for s=1:size(fir,2)
    fir(:,s)=(fir(:,s)-mean(fir(:,s)))/std(fir(:,s));
end

mskh=imread('avg_face_mask.tiff');
msk2=(mean(mskh,3)<127);

for k=1:800
    as(1:size(msk2,1),1:size(msk2,2))=0;
    for i=1:size(msk2,1)
        for j=1:size(msk2,2)
            if msk2(i,j)==1&&(i-asu(1,k))^2+(j-asu(2,k))^2<asu(3,k)^2
                as(i,j)=1;
            end
        end
    end
    asti(:,:,k)=as;
    astir(:,:,k)=as(:,end:-1:1);
end

for k=811:910
    as(1:size(msk2,1),1:size(msk2,2))=0;
    for i=1:size(msk2,1)
        for j=1:size(msk2,2)
            if msk2(i,j)==1&&(i-asu(1,k))^2+(j-asu(2,k))^2<asu(3,k)^2
                as(i,j)=1;
            end
        end
    end
    asti(:,:,k)=as;
    astir(:,:,k)=as(:,end:-1:1);
end

% compute fragment detection maps
for sz=1:size(fir,2)
    fira=fir(:,sz)';
    as=[];
    ap(1:size(msk2,1),1:size(msk2,2))=0;
    apl(1:size(msk2,1),1:size(msk2,2))=0;
    for k=1:800
        if leftright(sz)==0
            as=asti(:,:,k);
        else
            as=astir(:,:,k);
        end
        apl=apl+as;
        ap=ap+as*(fira(k)-mean(fira(801:810)));
    end
    ap=ap./apl;
    ap(~(ap<10000))=0;
    atmf(:,:,sz)=ap;
end

for k=1141:1940
    as(1:size(msk2,1),1:size(msk2,2))=0;
    for i=1:size(msk2,1)
        for j=1:size(msk2,2)
            if msk2(i,j)==1&&(i-asu(1,k))^2+(j-asu(2,k))^2<asu(3,k)^2
                as(i,j)=1;
            end
        end
    end
    asti(:,:,k)=as;
    astir(:,:,k)=as(:,end:-1:1);
end

% compute occlusion detection maps
for sz=1:size(fir,2)
    
    fira=fir(:,sz)';
    
    ap(1:size(msk2,1),1:size(msk2,2))=0;
    apl(1:size(msk2,1),1:size(msk2,2))=0;
    as=[];
    
    for k=1141:1940
        if leftright(sz)==0
            as=asti(:,:,k);
        else
            as=astir(:,:,k);
        end
        apl=apl+as;
        ap=ap+as*(fira(k)-mean(fira(911:1140)));
    end
    ap=ap./apl;
    ap(~(ap<10000))=0;
    avst2=ap;
    atmo(:,:,sz)=ap;
end

fsh=aps(1:3,911:1140);
fap=aps(4:6,911:1140);
aps(1:3,:)=aps(1:3,:)/sqrt(sum(std(fsh,0,2).^2))*sqrt(0.5);
aps(4:6,:)=aps(4:6,:)/sqrt(sum(std(fap,0,2).^2))*sqrt(0.5);

%compute fragment discrimination maps
for sz=1:size(fir,2)
    fira=fir(:,sz)';
    sta=aps(:,911:1140)*fira(911:1140)';
    generator=aps(:,:)'*sta;
    fr=fira;
    afr=[];
    numt=[];
    par=[];
    afr(1:size(asti,1),1:size(asti,2))=0;
    numt(1:size(asti,1),1:size(asti,2))=0.000000000000000000000000000001;
    par(1:size(asti,1),1:size(asti,2))=0;
    for i=1:1:800
        if leftright(sz)==0
            afr(1:size(asti,1),1:size(asti,2))=afr(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i)*fr(i);
            numt(1:size(asti,1),1:size(asti,2))=numt(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i);
            par(1:size(asti,1),1:size(asti,2))=par(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i)*generator(i);
        else
            afr(1:size(asti,1),1:size(asti,2))=afr(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i)*fr(i);
            numt(1:size(asti,1),1:size(asti,2))=numt(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i);
            par(1:size(asti,1),1:size(asti,2))=par(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i)*generator(i);
        end
    end
    meaf=afr./numt;
    meap=par./numt;
    asvt=[];
    bfr=[];
    bpar=[];
    avst(1:size(asti,1),1:size(asti,2))=0;
    bfr(1:size(asti,1),1:size(asti,2))=0;
    bpar(1:size(asti,1),1:size(asti,2))=0;
    for i=1:1:800
        if leftright(sz)==0
            bfr(1:size(asti,1),1:size(asti,2))=bfr(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf).*(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf);
            bpar(1:size(asti,1),1:size(asti,2))=bpar(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap).*(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap);
            avst(1:size(asti,1),1:size(asti,2))=avst(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf).*(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap);
        else
            bfr(1:size(asti,1),1:size(asti,2))=bfr(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf).*(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf);
            bpar(1:size(asti,1),1:size(asti,2))=bpar(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap).*(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap);
            avst(1:size(asti,1),1:size(asti,2))=avst(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf).*(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap);
        end
    end
    avst=avst./((bfr.*bpar).^0.5);
    avst(~(avst<1000))=0;
    avst1=avst;
    atm3(:,:,sz)=avst;
end


%compute occlusion discrimination maps

for sz=1:size(fir,2)
    fira=fir(:,sz)';
    sta=aps(:,911:1140)*fira(911:1140)';
    generator=aps(:,:)'*sta;
    fr=fira;
    afr=[];
    numt=[];
    par=[];
    afr(1:size(asti,1),1:size(asti,2))=0;
    numt(1:size(asti,1),1:size(asti,2))=0.000000000000000000000000000001;
    par(1:size(asti,1),1:size(asti,2))=0;
    for i=1141:1:1940
        if leftright(sz)==0
            afr(1:size(asti,1),1:size(asti,2))=afr(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i)*fr(i);
            numt(1:size(asti,1),1:size(asti,2))=numt(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i);
            par(1:size(asti,1),1:size(asti,2))=par(1:size(asti,1),1:size(asti,2))+asti(1:size(asti,1),1:size(asti,2),i)*generator(i);
        else
            afr(1:size(asti,1),1:size(asti,2))=afr(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i)*fr(i);
            numt(1:size(asti,1),1:size(asti,2))=numt(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i);
            par(1:size(asti,1),1:size(asti,2))=par(1:size(asti,1),1:size(asti,2))+astir(1:size(asti,1),1:size(asti,2),i)*generator(i);
        end
    end
    meaf=afr./numt;
    meap=par./numt;
    avst=[];
    bfr=[];
    bpar=[];
    avst(1:size(asti,1),1:size(asti,2))=0;
    bfr(1:size(asti,1),1:size(asti,2))=0;
    bpar(1:size(asti,1),1:size(asti,2))=0;
    for i=1141:1:1940
        if leftright(sz)==0
            bfr(1:size(asti,1),1:size(asti,2))=bfr(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf).*(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf);
            bpar(1:size(asti,1),1:size(asti,2))=bpar(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap).*(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap);
            avst(1:size(asti,1),1:size(asti,2))=avst(1:size(asti,1),1:size(asti,2))+(asti(1:size(asti,1),1:size(asti,2),i)*fr(i)-asti(1:size(asti,1),1:size(asti,2),i).*meaf).*(asti(1:size(asti,1),1:size(asti,2),i)*generator(i)-asti(1:size(asti,1),1:size(asti,2),i).*meap);
        else
            bfr(1:size(asti,1),1:size(asti,2))=bfr(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf).*(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf);
            bpar(1:size(asti,1),1:size(asti,2))=bpar(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap).*(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap);
            avst(1:size(asti,1),1:size(asti,2))=avst(1:size(asti,1),1:size(asti,2))+(astir(1:size(asti,1),1:size(asti,2),i)*fr(i)-astir(1:size(asti,1),1:size(asti,2),i).*meaf).*(astir(1:size(asti,1),1:size(asti,2),i)*generator(i)-astir(1:size(asti,1),1:size(asti,2),i).*meap);
            
        end
    end
    avst=avst./((bfr.*bpar).^0.5);
    avst(~(avst<1000))=0;
    avst1=avst;
    atm4(:,:,sz)=avst;
end

atmr=reshape(atmf,size(atmf,1)*size(atmf,2),size(atmf,3));
atmr(reshape((mean(astir,3).*mean(asti,3)==0),size(asti,1)*size(asti,2),1),:)=0;
atmf=reshape(atmr,size(atmf,1),size(atmf,2),size(atmf,3));
atmr=reshape(atmo,size(atmf,1)*size(atmf,2),size(atmf,3));
atmr(reshape((mean(astir,3).*mean(asti,3)==0),size(asti,1)*size(asti,2),1),:)=0;
atmo=reshape(atmr,size(atmf,1),size(atmf,2),size(atmf,3));
atmr=reshape(atm3,size(atmf,1)*size(atmf,2),size(atmf,3));
atmr(reshape((mean(astir,3).*mean(asti,3)==0),size(asti,1)*size(asti,2),1),:)=0;
atm3=reshape(atmr,size(atmf,1),size(atmf,2),size(atmf,3));
atmr=reshape(atm4,size(atmf,1)*size(atmf,2),size(atmf,3));
atmr(reshape((mean(astir,3).*mean(asti,3)==0),size(asti,1)*size(asti,2),1),:)=0;
atm4=reshape(atmr,size(atmf,1),size(atmf,2),size(atmf,3));

it4=wd_compit4_1;
figure;
subplot(2,4,1);
a1=mean(atmf,3)/max(max(abs(mean(atmf,3))))*0.5+0.5;
b1=[];
for i=1:size(a1,1)
    for j=1:size(a1,2)
        b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
    end
end
imshow(b1);
subplot(2,4,3);
a1=mean(atmo,3)/max(max(abs(mean(atmo,3))))*0.5+0.5;
b1=[];
for i=1:size(a1,1)
    for j=1:size(a1,2)
        b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
    end
end
imshow(b1);
subplot(2,4,5);
a1=mean(atm3,3)/max(max(abs(mean(atm3,3))))*0.5+0.5;
b1=[];
for i=1:size(a1,1)
    for j=1:size(a1,2)
        b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
    end
end
imshow(b1);
subplot(2,4,7);
a1=mean(atm4,3)/max(max(abs(mean(atm4,3))));
b1=[];
for i=1:size(a1,1)
    for j=1:size(a1,2)
        if a1(i,j)==0
            a1(i,j)=0.5;
        end
        b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
    end
end
imshow(b1);

aky=imread('average_face.tiff');
ap=mean(atmf,3);
ap(abs(ap)<=0.0001)=-10000;
apsi=sort(ap(:));
apsi=apsi(apsi>-1000);
a0=0.8;
subplot(2,4,2);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
ap=mean(atmo,3);
ap(abs(ap)<=0.0001)=10000;
apsi=sort(ap(:));
apsi=apsi(apsi<1000);
subplot(2,4,4);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
ap=mean(atm3,3);
ap(abs(ap)<=0.0001)=-10000;
apsi=sort(ap(:));
apsi=apsi(apsi>-1000);
subplot(2,4,6);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
ap=mean(atm4,3);
ap(abs(ap)<=0.0001)=10000;
apsi=sort(ap(:));
apsi=apsi(apsi<1000);
subplot(2,4,8);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);


atz=mean(abs(atmo),3);

for k=1:size(atmo,3)
    aty=atmo(:,:,k);
    atqmo(:,k)=aty(atz>0.00001);
    [I B]=sort(atqmo(:,k),'ascend');
    atqmoc(:,k)=(atqmo(:,k)<I(round(length(I)/5)));
    aty=atmf(:,:,k);
    atqmf(:,k)=aty(atz>0.00001);
    [I B]=sort(atqmf(:,k),'descend');
    atqmfc(:,k)=(atqmf(:,k)>I(round(length(I)/5)));
    aty=atm3(:,:,k);
    atqm3(:,k)=aty(atz>0.00001);
    [I B]=sort(atqm3(:,k),'descend');
    atqm3c(:,k)=(atqm3(:,k)>I(round(length(I)/5)));
    aty=atm4(:,:,k);
    atqm4(:,k)=aty(atz>0.00001);
end

overlap=mean(atqmoc&atqmfc,1)./mean(atqmoc|atqmfc,1);
overlap3c=mean(atqm3c&atqmfc,1)./mean(atqm3c|atqmfc,1);
overlap2=-diag(corr(atqmo,atqmf,'type','Spearman'));
overlap3=diag(corr(atqm3,atqmf,'type','Spearman'));

for i=1:3
    aky0(:,:,i)=double(aky)/255;
end
%superposition of two critical regions for consistent and inconsistent
%cells
overlapl=overlap;
[I B]=sort(overlapl,'descend');
a0=0.8;
figure;
for i=1:20
    akyy(1:size(aky,1),1:size(aky,2),1:3)=0;
    akyt=mean(akyy,3);
    ap=atmf(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*a0)))=1;
    akyy(:,:,1)=akyt;
    akyt(:,:)=0;
    ap=atmo(:,:,B(i));
    ap(abs(ap)<=0.0001)=10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi<1000);
    akyt(ap<apsi(round(length(apsi)*(1-a0))))=1;
    akyy(:,:,2)=akyt;
    subplot(4,10,i);
    aku=mean(akyy,3);
    for t=1:3
        akuu(:,:,t)=aku;
    end
    imshow(double(akyy)+0.6*double(aky0).*(akuu<0.0001));
end

for i=length(I)-19:length(I)
    akyy(1:size(aky,1),1:size(aky,2),1:3)=0;
    akyt=mean(akyy,3);
    ap=atmf(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*a0)))=1;
    akyy(:,:,1)=akyt;
    akyt(:,:)=0;
    ap=atmo(:,:,B(i));
    ap(abs(ap)<=0.0001)=10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi<1000);
    akyt(ap<apsi(round(length(apsi)*(1-a0))))=1;
    akyy(:,:,2)=akyt;
    subplot(4,10,i-(length(I)-20)+20);
    aku=mean(akyy,3);
    for t=1:3
        akuu(:,:,t)=aku;
    end
    imshow(double(akyy)+0.6*double(aky0).*(akuu<0.0001));
end


overlapl=overlap3c;
[I B]=sort(overlapl,'descend');
a0=0.8;
figure;
for i=1:20
    akyy(1:size(aky,1),1:size(aky,2),1:3)=0;
    akyt=mean(akyy,3);
    ap=atmf(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*a0)))=1;
    akyy(:,:,1)=akyt;
    akyt(:,:)=0;
    ap=atm3(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*(a0))))=1;
    akyy(:,:,2)=akyt;
    subplot(4,10,i);
    aku=mean(akyy,3);
    for t=1:3
        akuu(:,:,t)=aku;
    end
    imshow(double(akyy)+0.6*double(aky0).*(akuu<0.0001));
end

for i=length(I)-19:length(I)
    akyy(1:size(aky,1),1:size(aky,2),1:3)=0;
    akyt=mean(akyy,3);
    ap=atmf(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*a0)))=1;
    akyy(:,:,1)=akyt;
    akyt(:,:)=0;
    ap=atm3(:,:,B(i));
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    akyt(ap>apsi(round(length(apsi)*(a0))))=1;
    akyy(:,:,2)=akyt;
    subplot(4,10,i-(length(I)-20)+20);
    aku=mean(akyy,3);
    for t=1:3
        akuu(:,:,t)=aku;
    end
    imshow(double(akyy)+0.6*double(aky0).*(akuu<0.0001));
end

it4=wd_compit4_1;
overlap2p=overlap2;
[I B]=sort(overlap2p,'descend');
scol=[];
aky=imread('average_face.tiff');
scol(:,1)=mean(atqmf(:,B(1:40)),2);
scol(:,2)=mean(atqmo(:,B(1:40)),2);
figure;
for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    a1=aty/max(abs(aty(:)))*0.5+0.5;
    b1=[];
    for i=1:size(a1,1)
        for j=1:size(a1,2)
            b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
        end
    end
    subplot(1,8,2*k-1);
    imshow(b1);
end

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    ap=aty;
    
    if k~=2
        ap(abs(ap)<=0.0001)=-10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi>-1000);
        a0=0.8;
        subplot(1,8,2*k);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
    else
        ap(abs(ap)<=0.0001)=10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi<1000);
        subplot(1,8,2*k);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
    end
end

scol(:,1)=mean(atqmf(:,B(end-39:end)),2);
scol(:,2)=mean(atqmo(:,B(end-39:end)),2);

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    a1=aty/max(abs(aty(:)))*0.5+0.5;
    b1=[];
    for i=1:size(a1,1)
        for j=1:size(a1,2)
            b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
        end
    end
    subplot(1,8,2*k-1+4);
    imshow(b1);
end

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    ap=aty;
    
    if k~=2
        ap(abs(ap)<=0.0001)=-10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi>-1000);
        a0=0.8;
        subplot(1,8,2*k+4);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
    else
        ap(abs(ap)<=0.0001)=10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi<1000);
        subplot(1,8,2*k+4);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
    end
end


it4=wd_compit4_1;
ov3=overlap3;
[I B]=sort(ov3,'descend');
scol(:,1)=mean(atqmf(:,B(1:40)),2);
scol(:,2)=mean(atqm3(:,B(1:40)),2);
figure;
for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    a1=aty/max(abs(aty(:)))*0.5+0.5;
    b1=[];
    for i=1:size(a1,1)
        for j=1:size(a1,2)
            b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
        end
    end
    subplot(1,8,2*k-1);
    imshow(b1);
end

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    ap=aty;
    
    if k~=3
        ap(abs(ap)<=0.0001)=-10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi>-1000);
        a0=0.8;
        subplot(1,8,2*k);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
    else
        ap(abs(ap)<=0.0001)=10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi<1000);
        subplot(1,8,2*k);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
    end
end

scol(:,1)=mean(atqmf(:,B(end-39:end)),2);
scol(:,2)=mean(atqm3(:,B(end-39:end)),2);

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    a1=aty/max(abs(aty(:)))*0.5+0.5;
    b1=[];
    for i=1:size(a1,1)
        for j=1:size(a1,2)
            b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
        end
    end
    subplot(1,8,2*k-1+4);
    imshow(b1);
end

for k=1:2
    aty=atmo(:,:,1);
    aty(atz>0.00001)=scol(:,k);
    ap=aty;
    
    if k~=3
        ap(abs(ap)<=0.0001)=-10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi>-1000);
        a0=0.8;
        subplot(1,8,2*k+4);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
    else
        ap(abs(ap)<=0.0001)=10000;
        apsi=sort(ap(:));
        apsi=apsi(apsi<1000);
        subplot(1,8,2*k+4);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
    end
end
%tuning of consistent and inconsistent cells to fragment size 
facec=(mean(asti(:,:,811:910),3)>0);
for i=811:910
    facef=asti(:,:,i);
    fragr(i)=sum(facef(:)>0)/sum(facec(:));
end
[Ij Bj]=sort(fragr(811:910),'ascend');
xj(1)=0;
for p=1:5
    xj(p+1)=mean(Ij((p-1)*20+1:p*20));
end
xj(7)=1;
yj(1,:)=mean(fir(801:810,:),1);
for p=1:5
    yj(p+1,:)=mean(fir(Bj((p-1)*20+1:p*20)+810,:),1);
end
yj(7,:)=mean(fir(911:1140,:),1);
yj2=yj;
[I B]=sort(overlap2,'descend');
figure;
for j=1:2
    switch j
        case 1
            selb=B(1:40);
        case 2
            selb=B(end-39:end);
    end
    meat(:,j)=mean(yj2(:,selb),2);
    stdt(:,j)=std(yj2(:,selb),0,2)/sqrt(length(selb));
    stdt(:,j)=stdt(:,j)/(meat(7,j)-meat(1,j));
    meat(:,j)=(meat(:,j)-meat(1,j))/(meat(7,j)-meat(1,j));
    hold('all');
    errorbar(xj,meat(:,j),stdt(:,j));
end

%state space analysis
[coff2 scof2 latf2]=pca(fir(:,:));

for i=1:size(asti,1)
    for j=1:size(asti,2)
        xco(i,j)=j;
        yco(i,j)=i;
    end
end
xa=xco(atz>0.00001);
ya=yco(atz>0.00001);
for k=1:800
    xc0(k)=sum(sum(xco.*asti(:,:,k)))/sum(sum(asti(:,:,k)));
    yc0(k)=sum(sum(yco.*asti(:,:,k)))/sum(sum(asti(:,:,k)));
end
xc=(xc0-min(xc0))/(max(xc0)-min(xc0))-0.5;
yc=(yc0-min(yc0))/(max(yc0)-min(yc0))-0.5;
for k=1:800
    l0=[0.5+yc(k) 0.5-yc(k) 0.5];
    colo1(1:3,k)=l0;
    colo2(1:3,k)=[0.5+xc(k) 0.5-xc(k) 0.5];
    colo3(1:3,k)=[0.5+yc(k)-xc(k) 0.5-yc(k)-xc(k) 0.5+xc(k)];
    siz(k)=175;
end

xco2=(xco-min(xc0))/(max(xc0)-min(xc0))-0.5;
yco2=(yco-min(yc0))/(max(yc0)-min(yc0))-0.5;
for i=1:size(asti,1)
    for j=1:size(asti,2)
        if atz(i,j)>0.00001
            diagram(i,j,1:3)=[0.5+yco2(i,j)-xco2(i,j) 0.5-yco2(i,j)-xco2(i,j) 0.5+xco2(i,j)];
        else
            diagram(i,j,1:3)=1;
        end
    end
end
figure;imshow(diagram);



for i=1:size(aky,1)
    for j=1:size(aky,2)
        if ((i-167)/45)^2+((j-120)/85)^2<=1
            maskk(i,j)=1;
        else
            maskk(i,j)=0;
        end
    end
end


for i=1:size(atmf,3)
    ap=mean(atmf(:,:,i),3);
    ap(abs(ap)<=0.0001)=-10000;
    apsi=sort(ap(:));
    apsi=apsi(apsi>-1000);
    aq2(:,:,i)=(ap>apsi(round(length(apsi)*0.9)));
    asc(i)=sum(sum(aq2(:,:,i).*maskk==1))/sum(sum(aq2(:,:,i)));
end

[coff0, scof0 latf0]=pca(fir(1:800,asc>0.5));

figure;
subplot(2,1,1);
scatter(scof0(1:800,1),scof0(1:800,2),siz(1:800)*1.2,colo3(:,1:800)','Marker','.'); xlim([-9 16]);ylim([-9 11]);yticks([-10 -5 0 5 10]);
[coff0a scof0a latf0a]=pca(fir(1:800,:));
subplot(2,1,2);
scatter(scof0a(1:800,1),scof0a(1:800,2),siz(1:800)*1.2,colo3(:,1:800)','Marker','.'); xlim([-9 16]);ylim([-9 11]);yticks([-10 -5 0 5 10]);

scof3=scof2;
[I B]=sort(fragr(811:910));
size_step=[];
for k=1:5
    size_step(k+1,:)=mean(scof3(810+B((k-1)*20+1:k*20),:),1);
end
size_step(1,:)=mean(scof3(801:810,:),1);
size_step(size(size_step,1)+1,:)=mean(scof3(911:1200-60,:),1);
figure;
scatter(scof3(1:800,1),scof3(1:800,3),'MarkerFaceColor',[1 0.83 0.83],'MarkerEdgeColor',[1 0.83 0.83]);hold('all');
scatter(scof3(1141:1940,1),scof3(1141:1940,3),'MarkerFaceColor',[0.83 0.83 1],'MarkerEdgeColor',[0.83 0.83 1]);
scatter(scof3(801:1140,1),scof3(801:1140,3),'MarkerFaceColor',[0.6 0.6 0.6],'MarkerEdgeColor',[0.6 0.6 0.6]);
plot(size_step(:,1),size_step(:,3),'color',[0 0 0]','linewidth',2);
scatter(size_step(:,1),size_step(:,3),'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);

figure;
scatter(scof3(1:800,1),scof3(1:800,2),'MarkerFaceColor',[1 0.83 0.83],'MarkerEdgeColor',[1 0.83 0.83]);hold('all');
scatter(scof3(1141:1940,1),scof3(1141:1940,2),'MarkerFaceColor',[0.83 0.83 1],'MarkerEdgeColor',[0.83 0.83 1]);
scatter(scof3(801:1140,1),scof3(801:1140,2),'MarkerFaceColor',[0.6 0.6 0.6],'MarkerEdgeColor',[0.6 0.6 0.6]);
plot(size_step(:,1),size_step(:,2),'color',[0 0 0]','linewidth',2);
scatter(size_step(:,1),size_step(:,2),'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);


for theta=0:0.01:2*pi+0.01
    scof4=scof3(:,1)*cos(theta)+scof3(:,3)*sin(theta);
    dp(round(theta*100+1))=(mean(scof4(1:800,1))-mean(scof4(911:1140,1)))/sqrt(((std(scof4(1:800,1)).^2)*799+(std(scof4(911:1140,1)).^2)*229)/(799+229));
    xii(round(theta*100+1))=abs(dp(round(theta*100+1)))*cos(theta);
    yii(round(theta*100+1))=abs(dp(round(theta*100+1)))*sin(theta);
    lz10(round(theta*100+1))=std(scof4(1:800,1));
    lz20(round(theta*100+1))=std(scof4(911:1140,1));    
end
figure;polarplot(0:0.01:2*pi+0.01,abs(dp));
astirs=reshape(asti(:,:,:),size(asti,1)*size(asti,2),size(asti,3));
sco3(:,1)=scof3(:,1);
size_step=[];
projd=[];
for p=2:50
    projd(1:100,1)=0;
    projd(p)=1;
    pl(1:100,p)=projd/norm(projd);
    sco3(:,2)=scof3(:,2:100)*projd(2:end)/norm(projd(2:end));
    dp=[];
    for theta=0:pi/360:pi/360*359
        firp=sco3(:,1)*cos(theta)+sco3(:,2)*sin(theta);
        dp(round(theta/(pi/360)+1))=(mean(firp(1:800,1))-mean(firp(911:1140,1)))/sqrt(((std(firp(1:800,1)).^2)*799+((std(firp(911:1140,1))).^2)*229)/(799+229));
    end
    atmf2=astirs(:,1:800)*(sco3(1:800,2)-mean(sco3(801:810,2)));
    atmo2=astirs(:,1141:1940)*(sco3(1141:1940,2)-mean(sco3(911:1140,2)));
    cosmp(p)=-corr(atmf2(atz(:)>0.00001),atmo2(atz(:)>0.00001));
    dp0p(p)=max(abs(dp));
    size_step=[];
    size_step(1,:)=mean(sco3(801:810,:),1);
    for k=1:5
     size_step(k+1,:)=mean(sco3(810+B((k-1)*20+1:k*20),:),1);
    end
    size_step(7,:)=mean(sco3(911:1200-60,:),1);
    lin_idp(p)=corr(size_step(:,2),xj');
    sis(:,1)=size_step(:,1)-size_step(7,1);
    sis(:,2)=size_step(:,2)-size_step(7,2);
    curvatp(p)=abs(sis(1,1)*sis(2,2)-sis(2,1)*sis(1,2)-sis(3,1)*sis(2,2)+sis(2,1)*sis(3,2)+sis(3,1)*sis(4,2)-sis(4,1)*sis(3,2)-sis(5,1)*sis(4,2)+sis(4,1)*sis(5,2)+sis(5,1)*sis(6,2)-sis(6,1)*sis(5,2))/2/(norm(sis(1,:))^2);
end

figure;
subplot(2,1,1);
plot(2:50,curvatp(2:50));
subplot(2,1,2);
plot(2:50,dp0p(2:50));

figure;
for kn=1:3
as=[];
ap=[];
apl=[];
ap(1:size(msk2,1),1:size(msk2,2))=0;
apl(1:size(msk2,1),1:size(msk2,2))=0;
fira=scof3(:,kn);     
for k=1:800
    as=asti(:,:,k);
    apl=apl+as;
    ap=ap+as*(fira(k)-mean(fira(801:810)));
end
apt(:,:,1)=ap./apl;
ap(1:size(msk2,1),1:size(msk2,2))=0;
apl(1:size(msk2,1),1:size(msk2,2))=0;
as=[];
for k=1141:1940
    as=asti(:,:,k);
    apl=apl+as;
    ap=ap+as*(fira(k)-mean(fira(911:1140)));
end
apt(:,:,2)=ap./apl;
apu=apt(:,:,1);
apv=apt(:,:,2);
cop(kn)=corr(apu(apu<1000),apv(apu<1000),'Type','Spearman');
for k=1:2
   aty=apt(:,:,k);
   a1=aty/max(abs(aty(:)))*0.5+0.5;
   a1(~(a1<1000))=0.5;
   b1=[];
   for i=1:size(a1,1)
    for j=1:size(a1,2)
      b1(i,j,:)=it4(round(a1(i,j)*255)+1,:);
    end
   end
   subplot(3,4,2*k+3+(kn-2)*4);
   imshow(b1);
end
for k=1:2
   aty=apt(:,:,k);
   aty(~(aty<10000))=0;
   ap=aty;
   
   if k==1
   ap(abs(ap)<=0.0001)=-10000;
   apsi=sort(ap(:));
   apsi=apsi(apsi>-1000);
   a0=0.8;
   subplot(3,4,2*k+4+(kn-2)*4);imshow(0.2*double(aky)+0.8*double(aky).*(ap>apsi(round(length(apsi)*a0))),[]);
   else
   ap(abs(ap)<=0.0001)=10000;
   apsi=sort(ap(:));
   apsi=apsi(apsi<1000);
   subplot(3,4,2*k+4+(kn-2)*4);imshow(0.2*double(aky)+0.8*double(aky).*(ap<apsi(round(length(apsi)*(1-a0)))),[]);
   end
end
end

%population decoding analysis
colorcode=[0 0 1;0 1 1; 0 0.5 0; 1 0 0; 1 0.7 0; 1 0 0.7];
figure;
featc=[];
for s=1:6
for i=1:30
    g2=scof2(911:1140,1:i);
    g2(:,end+1)=1;
    for t=1:10
        selecf=[1:23*(t-1) 23*t+1:230];
        [b bint r rint stats]=regress(aps(s,910+selecf)',g2(selecf,:));
        complement=[23*(t-1)+1:23*t];
        pred(complement)=g2(complement,:)*b;
    end
    l1=pred(1:230)';
    l2=aps(s,911:1140)';
    evf(s,i)=1-sum((l1-l2).^2)/sum((l2-mean(l2)).^2);
    corf(s,i)=corr(l1,l2,'type','Spearman');
    featc(:,s,i)=l1;
    if s==4
    subplot(1,4,1);
    scatter(l1,l2);
    end
end
end
for i=1:30
   a1=aps(:,911:1140)';
   b1=squeeze(featc(:,:,i));
   for k=1:size(a1,1)
       a1(k,:)=a1(k,:)/norm(a1(k,:));
       b1(k,:)=b1(k,:)/norm(b1(k,:));
   end
   cosa=a1*b1';
   apn=[];
   for k=1:size(a1,1)
       apn(k)=sum(cosa(:,k)>=cosa(k,k));
   end
   erate_com(i)=(mean(apn)-1)/(size(a1,1)-1);
end
subplot(1,4,2);
for i=1:6
    hold('all');
plot(corf(i,1:30)','MarkerEdgeColor',colorcode(i,:),'Color',colorcode(i,:),'MarkerFaceColor',[1 1 1],'Marker','o');
end
xlim([-1 32]);

featc=[];
pred=[];
for s=1:6
for i=1:50
    g2=scof2(1:800,1:i);
    g2(:,end+1)=1;
    for t=1:10
        selecf=[1:80*(t-1) 80*t+1:800];
        [b bint r rint stats]=regress(aps(s,selecf)',g2(selecf,:));
        complement=[80*(t-1)+1:80*t];
        pred(complement)=g2(complement,:)*b;
    end
    l1=pred(1:800)';
    l2=aps(s,1:800)';
    featc(:,s,i)=l1;
    evff(s,i)=1-sum((l1-l2).^2)/sum((l2-mean(l2)).^2);
    corff(s,i)=corr(l1,l2,'type','Spearman');   
end
end
for i=1:30
   a1=aps(:,1:800)';
   b1=squeeze(featc(:,:,i));
   for k=1:size(a1,1)
       a1(k,:)=a1(k,:)/norm(a1(k,:));
       b1(k,:)=b1(k,:)/norm(b1(k,:));
   end
   cosa=a1*b1';
   apn=[];
   for k=1:size(a1,1)
       apn(k)=sum(cosa(:,k)>=cosa(k,k));
   end
   erate_fra(i)=(mean(apn)-1)/(size(a1,1)-1);
end

subplot(1,4,4);
for i=1:6
    hold('all');
plot(corff(i,1:30)','MarkerEdgeColor',colorcode(i,:),'Color',colorcode(i,:),'MarkerFaceColor',[1 1 1],'Marker','o');
end
xlim([-1 32]);

featc=[];
pred=[];
for s=1:6
for i=1:30
    g2=scof2(1141:1940,1:i);
    g2(:,end+1)=1;
    for t=1:10
        selecf=[1:80*(t-1) 80*t+1:800];
        [b bint r rint stats]=regress(aps(s,selecf+1140)',g2(selecf,:));
        complement=[80*(t-1)+1:80*t];
        pred(complement)=g2(complement,:)*b;
    end
    l1=pred(1:800)';
    l2=aps(s,1141:1940)';
    featc(:,s,i)=l1;    
    evof(s,i)=1-sum((l1-l2).^2)/sum((l2-mean(l2)).^2);
    corof(s,i)=corr(l1,l2,'type','Spearman');
end
end

for i=1:30
   a1=aps(:,1141:1940)';
   b1=squeeze(featc(:,:,i));
   for k=1:size(a1,1)
       a1(k,:)=a1(k,:)/norm(a1(k,:));
       b1(k,:)=b1(k,:)/norm(b1(k,:));
   end
   cosa=a1*b1';
   apn=[];
   for k=1:size(a1,1)
       apn(k)=sum(cosa(:,k)>=cosa(k,k));
   end
   erate_ocl(i)=(mean(apn)-1)/(size(a1,1)-1);
end

subplot(1,4,3);
for i=1:6
    hold('all');
plot(corof(i,1:30)','MarkerEdgeColor',colorcode(i,:),'Color',colorcode(i,:),'MarkerFaceColor',[1 1 1],'Marker','o');
end
xlim([-1 32]);


colorcode=[0 0 1;0 1 1; 0 0.5 0; 1 0 0; 1 0.7 0; 1 0 0.7];
featc=[];
figure;

featc=[];
for s=1:6
for i=1:30
    g2=scof2(1141:1940,1:i);
    g2(:,end+1)=1;
    g02=scof2(911:1140,1:i);
    g02(:,end+1)=1;    
    [b bint r rint stats]=regress(aps(s,1141:1940)',g2(1:800,:));
    l1=(g02*b);
    l2=aps(s,911:1140)';
    featc(:,s,i)=l1;       
    evo2c(s,i)=1-sum((l1-l2).^2)/sum((l2-mean(l2)).^2);
    coro2c(s,i)=corr(l1,l2,'type','Spearman'); 
    if s==4
     subplot(1,4,1);
     scatter(l1,l2,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
    end
end
end

for i=1:30
   a1=aps(:,911:1140)';
   b1=squeeze(featc(:,:,i));
   for k=1:size(a1,1)
       a1(k,:)=a1(k,:)/norm(a1(k,:));
       b1(k,:)=b1(k,:)/norm(b1(k,:));
   end
   cosa=a1*b1';
   apn=[];
   for k=1:size(a1,1)
       apn(k)=sum(cosa(:,k)>=cosa(k,k));
   end
   erate_ocl2com(i)=(mean(apn)-1)/(size(a1,1)-1);
end

subplot(1,4,2);
for i=1:6
    hold('all');
plot(coro2c(i,1:30)','MarkerEdgeColor',colorcode(i,:),'Color',colorcode(i,:),'MarkerFaceColor',[1 1 1],'Marker','o');
end
xlim([-2 32]);

featc=[];
for s=1:6
for i=1:30
    g2=scof2(1:800,1:i);
    g2(:,end+1)=1;
    g02=scof2(911:1140,1:i);
    g02(:,end+1)=1;    
    [b bint r rint stats]=regress(aps(s,1:800)',g2(1:800,:));
    l1=(g02*b);
    l2=aps(s,911:1140)';
    featc(:,s,i)=l1;      
    evf2c(s,i)=1-sum((l1-l2).^2)/sum((l2-mean(l2)).^2);
    corf2c(s,i)=corr(l1,l2,'type','Spearman'); 
    if s==4
     subplot(1,4,3);
    scatter(l1,l2,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
    end
end
end


subplot(1,4,4);
for i=1:6
    hold('all');
plot(corf2c(i,1:30)','MarkerEdgeColor',colorcode(i,:),'Color',colorcode(i,:),'MarkerFaceColor',[1 1 1],'Marker','o');
end
xlim([-2 32]);

for i=1:30
   a1=aps(:,911:1140)';
   b1=squeeze(featc(:,:,i));
   for k=1:size(a1,1)
       a1(k,:)=a1(k,:)/norm(a1(k,:));
       b1(k,:)=b1(k,:)/norm(b1(k,:));
   end
   cosa=a1*b1';
   apn=[];
   for k=1:size(a1,1)
       apn(k)=sum(cosa(:,k)>=cosa(k,k));
   end
   erate_fra2com(i)=(mean(apn)-1)/(size(a1,1)-1);
end
figure;
plot(1-erate_com);
hold('all');
plot(1-erate_fra2com);
plot(1-erate_ocl2com);

%population representation of "identity", face completeness and condition
newc=aps([1 4],:)';
t=0;
tos(1:2,1:2)=[1 2;4 3];
ls=[];
lt1=[];lt2=[];
for i=1:2
    for j=1:2
        t=t+1;
        lt1(1:800)=0;
        if i==1
        lt1(newc(1:800,1)<=0)=1;
        else
        lt1(newc(1:800,1)>0)=1;
        end
        lt2(1:800)=0;
        if j==1
          lt2(newc(1:800,2)<=0)=1; 
        else
          lt2(newc(1:800,2)>0)=1; 
        end
        fir2o=fir(1:800,:);
        ls(:,tos(i,j))=mean(fir(lt1==1&lt2==1,:),1);
    end
end
lt1=[];lt2=[];
for i=1:2
    for j=1:2
        t=t+1;
        lt1(1:800)=0;
        if i==1
        lt1(newc(1141:1940,1)<=0)=1;
        else
        lt1(newc(1141:1940,1)>0)=1;
        end
        lt2(1:800)=0;
        if j==1
          lt2(newc([1141:1940],2)<=0)=1; 
        else
          lt2(newc([1141:1940],2)>0)=1; 
        end 
        fir2o=fir(1141:1940,:);
        lso(:,tos(i,j))=mean(fir2o(lt1==1&lt2==1,:),1);         
    end
end
lt1=[];lt2=[];
for i=1:2
    for j=1:2
        t=t+1;
        lt1(1:230)=0;
        if i==1
        lt1(newc(911:1140,1)<=0)=1;
        else
        lt1(newc(911:1140,1)>0)=1;
        end
        lt2(1:230)=0;
        if j==1
          lt2(newc(911:1140,2)<=0)=1; 
        else
          lt2(newc(911:1140,2)>0)=1; 
        end  
        fir2o=fir(911:1140,:);
        lsc(:,tos(i,j))=mean(fir2o(lt1==1&lt2==1,:),1);           
    end
end

[coffr scofr latfr]=pca(ls');
scofr=ls'*coffr;
scofro=lso'*coffr;
scofrc=lsc'*coffr;
for i=1:4
    sizecode(i)=500;
end
colorcode=[1 0 0; 1 0.5 0; 0 1 0; 0 0 1];
figure;
subplot(1,3,1);hold('all');
plot(scofro([1:end 1],1),scofro([1:end 1],2),'color',[0.6 0.6 0.6],'linewidth',1.2);plot(scofrc([1:end 1],1),scofrc([1:end 1],2),'color',[0 0 0],'linewidth',1.2);
scatter(scofro(:,1),scofro(:,2),sizecode,colorcode,'Marker','.');hold('all');
scatter(scofrc(:,1),scofrc(:,2),sizecode,colorcode,'Marker','.');hold('all');plot(scofr([1:end 1],1),scofr([1:end 1],2),'color',[0 0 0],'linewidth',1.2,'LineStyle','--');scatter(scofr(:,1),scofr(:,2),sizecode,colorcode,'Marker','.');

[coffr scofr latfr]=pca(lso');
scofr=ls'*coffr;
scofro=lso'*coffr;
scofrc=lsc'*coffr;
subplot(1,3,2);hold('all');
plot(scofro([1:end 1],1),scofro([1:end 1],2),'color',[0.6 0.6 0.6],'linewidth',1.2);plot(scofrc([1:end 1],1),scofrc([1:end 1],2),'color',[0 0 0],'linewidth',1.2);
scatter(scofro(:,1),scofro(:,2),sizecode,colorcode,'Marker','.');hold('all');
scatter(scofrc(:,1),scofrc(:,2),sizecode,colorcode,'Marker','.');hold('all');plot(scofr([1:end 1],1),scofr([1:end 1],2),'color',[0 0 0],'linewidth',1.2,'LineStyle','--');scatter(scofr(:,1),scofr(:,2),sizecode,colorcode,'Marker','.');

[coffr scofr latfr]=pca(lsc');
scofr=ls'*coffr;
scofro=lso'*coffr;
scofrc=lsc'*coffr;
subplot(1,3,3);hold('all');plot(scofro([1:end 1],1),scofro([1:end 1],2),'color',[0.6 0.6 0.6],'linewidth',1.2);plot(scofrc([1:end 1],1),scofrc([1:end 1],2),'color',[0 0 0],'linewidth',1.2);
scatter(scofro(:,1),scofro(:,2),sizecode,colorcode,'Marker','.');hold('all');
scatter(scofrc(:,1),scofrc(:,2),sizecode,colorcode,'Marker','.');hold('all');plot(scofr([1:end 1],1),scofr([1:end 1],2),'color',[0 0 0],'linewidth',1.2,'LineStyle','--');scatter(scofr(:,1),scofr(:,2),sizecode,colorcode,'Marker','.');

lscon(:,1)=mean(fir(801:810,:),1);
lscon(:,2)=mean(fir(1:800,:),1);
lscon(:,4)=mean(fir(911:1140,:),1);
lscon(:,3)=mean(fir(1141:1940,:),1);
[coffrc scofrc latfrc]=pca(lscon');

[Ij Bj]=sort(fragr(811:910),'ascend');
lscomp(:,1)=mean(fir(801:810,:),1);
for i=1:5
    lscomp(:,i+1)=mean(fir(Bj((i-1)*20+1:i*20)+810,:),1);
end
lscomp(:,7)=mean(fir(911:1140,:),1);
[coffrcom scofrcom latfrcom]=pca(lscomp');

figure;
subplot(1,2,1);scatter(scofrc(:,1),scofrc(:,2));hold('all');plot(scofrc(:,1),scofrc(:,2));
subplot(1,2,2);scatter(scofrcom(:,1),scofrcom(:,2));hold('all');plot(scofrcom(:,1),scofrcom(:,2));
