function gamutMappingStep3

%gamutMappingStep3 uses soft compression to the gamut center with lambda=2/3 to compress 
%any source pixels that are outside of destination gamut into destination 
%gamut which is done by 5 degree hue sector. 
%In each hue sector, center compression is done by 10 degree sub hue
%sector and center is at 1/2 of destination lightness
%Also it allows you to display 2D hue sectors view after center compression.

    close all
    load('data/keyInfo.mat')

    index=1;
    %Gamut Mapping
    for j=-180:5:175


        %Code below is the same from step 1 & step 2
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        i1=find(hue_dest>j);%find hue angle larger than j in destination gamut
        i2=find(hue_dest(i1)<=j+5);%find hue angle equal or smaller than j+5 in destination gamue
        dest_list=i1(i2);%dest is list of Indigo/destination gamut point inside this sector
        ll=length(dest_list);
        dest_list(ll+1:ll+2)=[1 num1];%add white and black from destination gamut to this list

        i1=find(hue_src>j);%find hue angle larger than j in source gamut
        i2=find(hue_src(i1)<=j+5);%find hue angle equal or smaller than j+5 in source gamut
        src_list=i1(i2);%src is list of Image/source gamut pixels inside this sector

        %for Destination gamut
        ch_dest_list=ch_dest(dest_list);%Chroma values in this sector
        ch_dest_max=max(ch_dest(dest_list));%max chroma value in this sector

        l_dest_list=l_dest(dest_list);%lightness values in this sector
        l_dest_max=max(l_dest(dest_list));%max chroma value in this sector
        l_dest_min=min(l_dest(dest_list));%min chroma value in this sector



        x=[0,0,ch_dest_max,ch_dest_max,0];
        y=[l_dest_max,l_dest_min,l_dest_min,l_dest_max,l_dest_max];


        %for Image/Source gamut
        [nrow,~]=size(src_list);
        ch_src_list=ch_src(src_list);%Chroma values in this sector
        ch_src_max=max(ch_src(src_list));% max chroma value in this sector

        l_src_list=l_src(src_list);%all lightness values in this sector
        l_src_max=max(l_src(src_list));%max lightness value in this sector

        %Compress Chroma into Bounding Cylinder
        lambda = 2/3;
        if ch_src_max > ch_dest_max
            for i=1:nrow
                if ch_src_list(i) > ch_dest_max
                    ch_src_list(i) =  (1-lambda)*(ch_dest_max/ch_src_max)*ch_src_list(i) + lambda*ch_dest_max;
                else
                    ch_src_list(i) =  (1-lambda)*(ch_dest_max/ch_src_max)*ch_src_list(i) + lambda*ch_src_list(i);
                end

            end
        end
        %use MATLAB build-in function 'convhull' to form a gamut convex hull from
        %the chroma and lightness for this sector
        dest_list_index = convhull(ch_dest(dest_list),l_dest(dest_list));
%         src_list_index = convhull(ch_src_list,l_src(src_list));
        if length(src_list) > 5
            src_list_index = convhull(ch_src(src_list),l_src(src_list));
        else
            src_list_index = [1:length(src_list)];
        end

        l_img_list_ref=l_src_list;



        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Code above is the same from step 1 & step 2




        %Code below is step 3: Center compression

        %Prepare For Center compression
        l_dest_list=l_dest_list-center;%Shift down so that mid point moved to origin
        l_dest_list_max=max(l_dest_list);%max lightness value in destination gamut
        l_dest_list_min=min(l_dest_list);%min lightness value in destination gamut
        l_d_convex=l_dest_list(dest_list_index); %convex list in this hue sector
        [l_d_convex_unique,ia2,~]=unique(l_d_convex);
        ch_d_convex=ch_dest_list(dest_list_index);
        lenth_l_d_convex_unique=length(l_d_convex_unique);
        l_src_list=l_src_list-center;% source lightness values in this hue sector
        l_img_list_ref=l_img_list_ref-center;%save another copy



        %Calculate sector_lightness and sector_chroma for destination gamut
        %in this hue sector
        %In each hue sector, horizontally divide the destination gamut into length of 0.05
        %sections and find the maximum chroma value in each gamut
        %section. The maximum chroma is what we can achieve when printed
        %This step is the preparion of center compression
        index2=1;
        for light=l_dest_list_min:0.05:l_dest_list_max
            sector_lightness(index2)=light;
            index2=index2+1;
        end
        sector_lightness(index2)=l_dest_list_max;
        for i=1:index2
            for m=1:lenth_l_d_convex_unique-1
                if (sector_lightness(i) >= l_d_convex_unique(m)) && (sector_lightness(i) <= l_d_convex_unique(m+1))
                    sector_chroma(i)=(sector_lightness(i)-l_d_convex_unique(m))*(ch_d_convex(ia2(m+1))-ch_d_convex(ia2(m)))/(l_d_convex_unique(m+1)-l_d_convex_unique(m))+ch_d_convex(ia2(m));
                end
            end
        end


        sub_hue_sector_dest=zeros(index2,1);
        sat_inv_sector_dest=zeros(index2,1);
        %%%%%%%%%%%%%%%%%%%%%%%
        for i = 1:index2
            %sub_hue_sector_dest is angle between destination lightness and destination chroma of a given
            %pixel in destination gamut
            hue = (atan(abs(sector_lightness(i))/abs(sector_chroma(i)))/pi)*180;
            if sector_lightness(i)==0
                sub_hue_sector_dest(i) = 0;
            elseif sector_lightness(i)>=0 % First
                sub_hue_sector_dest(i) = hue;
            elseif sector_lightness(i)<=0 % Fourth
                sub_hue_sector_dest(i) = - hue;
            end
            %sat_inv_sector_dest = sqrt((destination lightness)^2 + (destination chroma)^2)
            %which is the Euclidean distance from origin of a pair of lightness and chroma
            sat_inv_sector_dest(i)=sqrt((sector_lightness(i))^2+(sector_chroma(i))^2);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        sub_hue_sector_src=zeros(index2,1);
        sat_inv_sector_src=zeros(index2,1);
        sat_inv_sector_src_ref=zeros(index2,1);
        %%%%Calculate sub_hue_sector and saturation inverse for source gamut
        for i = 1:nrow
            %hue_sector is angle between lightness and chroma of a given
            %pixel in source gamut
            hue = (atan(abs(l_src_list(i))/abs(ch_src_list(i)))/pi)*180;
            if l_src_list(i)==0
                sub_hue_sector_src(i) = 0;
            elseif l_src_list(i)>=0 % First
                sub_hue_sector_src(i) = hue;
            elseif l_src_list(i)<=0 % Fourth
                sub_hue_sector_src(i) = - hue;
            end
            %sat_inv_sector_src = sqrt((source lightness)^2 + (source chroma)^2)
            %which is the Euclidean distance from origin of a pair of lightness and chroma
            sat_inv_sector_src(i)=sqrt((l_src_list(i))^2+(ch_src_list(i))^2); 
            sat_inv_sector_src_ref(i)=sqrt((l_src_list(i))^2+(ch_src_list(i))^2);% save another copy
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%for every sub Hue_sector angle
        for jindex=-90:10:80
            index1=find(sub_hue_sector_dest>jindex);%find all pixels in destination gamut where sub_hue_sector angle is greater than jindex
            index2=find(sub_hue_sector_dest(index1)<=jindex+10);% find all pixels in destination gamut where sub_hue_sector angle is less than jindex+10
            A2=index1(index2); %list of pixels of this sub hue sector of destination

            sat_inv_sector_dest_max=min(sat_inv_sector_dest(A2));

            index3=find(sub_hue_sector_src>jindex);%find all pixels in source gamut where sub_hue_sector angle is greater than jindex
            index4=find(sub_hue_sector_src(index3)<=jindex+10);% find all pixels in source gamut where sub_hue_sector angle is less than jindex+10
            A3=index3(index4);%list of pixels of this sub hue sector of source
            length_A3=length(A3);
            sat_inv_sector_src_list=sat_inv_sector_src(A3);
            sat_inv_sector_src_max=max(sat_inv_sector_src(A3));

            lambda=2/3;
            % Center Compression of sat_inv_sector_src
            if sat_inv_sector_src_max > sat_inv_sector_dest_max
                for i=1:length_A3
                    if sat_inv_sector_src_list(i) > sat_inv_sector_dest_max
                        sat_inv_sector_src_list(i) =  (1-lambda)*(sat_inv_sector_dest_max/sat_inv_sector_src_max)*sat_inv_sector_src_list(i) + lambda*sat_inv_sector_dest_max;
                    else
                        sat_inv_sector_src_list(i) =  (1-lambda)*(sat_inv_sector_dest_max/sat_inv_sector_src_max)*sat_inv_sector_src_list(i) + lambda*sat_inv_sector_src_list(i);
                    end

                end
            end
            sat_inv_sector_src(A3)=sat_inv_sector_src_list; %%%Update calculation
        end
        
        %Map sat_inv_sector_src back to lightness and chroma
        for i = 1:nrow
            scale = 0.97 * sat_inv_sector_src(i)/sat_inv_sector_src_ref(i);
            l_src_list(i) = l_src_list(i) * scale;
            ch_src_list(i) = ch_src_list(i) * scale;

        end

        l_src_list = l_src_list + center;%shift source gamut up by 1/2 of destination lightness
        ch_src(src_list)=ch_src_list; %Update changes to original chroma list
        l_src(src_list)=l_src_list; %Update changes to original lightness list

        %use MATLAB build-in function 'convhull' to form a gamut convex hull from
        %the chroma and lightness for this sector
        dest_list_index = convhull(ch_dest(dest_list),l_dest(dest_list));
        if length(src_list) > 5
            src_list_index = convhull(ch_src(src_list),l_src(src_list));
        else
            src_list_index = [1:length(src_list)];
        end
        
        

        %Comment out plotting to speed process,
        %this plot is showing Lightness v.s. Chroma after center compression

        %plot
        figure
        %%%Plot limit
        plot(x,y)
        hold on
        %%%Indigo
        plot(ch_dest(dest_list),l_dest(dest_list),'.','MarkerSize',10)
        hold on
        plot(ch_dest(dest_list(dest_list_index)),l_dest(dest_list(dest_list_index)),'o-','MarkerSize',10)
        hold on
        %%%Image
        plot(ch_src_list,l_src_list,'.','MarkerSize',10)
        hold on
        plot(ch_src_list(src_list_index),l_src_list(src_list_index),'o-','MarkerSize',10)

        xlim([0 180])
        ylim([0 220])
        grid on
        legend('Bounding Cylinder','pixels of Destination Gamut in this hue sector','Convex hull of pixels of Destination Gamut','pixels of Source Gamut in this hue sector','Convex hull of pixels of Source Gamut')
        xlabel('Chroma');ylabel('Lightness')
        input = sprintf('Step 3 After center compression slice view (%d,%d]',j,j+5);
        title(input)
        name = sprintf('Step 3 After center compression slice view (%d,%d].png',j,j+5);
        saveas(gca,[pwd '/output/hueSlicePlots/',name]);
        close all


        index=index+1;

    end


    save data/ch_src_gamutMapped.mat ch_src
    save data/l_src_gamutMapped.mat l_src

end



