S1_ExploitGroupSparcityToGenerateOutputImage.m

%   Copyright 2010 Chih-Yuan Yang, Jia-Bin Huang, and Ming-Hsuan Yang
%
%   Licensed under the Apache License, Version 2.0 (the "License");
%   you may not use this file except in compliance with the License.
%   You may obtain a copy of the License at
%
%       http://www.apache.org/licenses/LICENSE-2.0
%
%   Unless required by applicable law or agreed to in writing, software
%   distributed under the License is distributed on an "AS IS" BASIS,
%   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%   See the License for the specific language governing permissions and
%   limitations under the License.
%Chih-Yuan Yang, EECS, UC Merced
%2/22/2013

Para.bExtractFeature = true;    %I think that it may be better to hide this two debugging variables. They should be true in most cases.
Para.bDoCluster = true;
LayerNumber = length(AllLayers);

if Para.Zooming == 3
    PairNumber = 6;
    LowLayerIdx_Array = [11 10 9 8 7 6];
    HighLayerIdx_Array = [6 5 4 3 2 1];
    LayerIdx_I0 = 6;
elseif Para.Zooming == 4
    PairNumber = 7;
    LowLayerIdx_Array = [14 13 12 11 10 9 8 ];
    HighLayerIdx_Array = [7 6 5 4 3 2 1];
    LayerIdx_I0 = 8;
else
    error('scaling factor not supported');
end
img_input_y = AllLayers(LayerIdx_I0).GridAsHighLayer;

if Para.bExtractFeature
    %build the rebuild-high
    MappingLayerLow(PairNumber,1) = struct('Index', 0 , 'TrueWidth', 0 , 'TrueHeight' , 0 , 'FormatWidth' , 0 , 'FormatHeight' , 0 , ...
        'ValidWidth' , 0 , 'ValidHeight' , 0, 'Image' , []);
    MappingLayerHigh(PairNumber,1) = struct('Index', 0 , 'TrueWidth', 0 , 'TrueHeight' , 0 , 'FormatWidth' , 0 , 'FormatHeight' , 0 , ...
        'ValidWidth' , 0 , 'ValidHeight' , 0, 'Image' , []);

    for i=1:PairNumber           
        LowLayerIdx = LowLayerIdx_Array(i);
        LowLayer = AllLayers(LowLayerIdx);          %here may not be 14, depends on the Zooming parameter.
        HighLayerIdx = HighLayerIdx_Array(i);
        OverBigLayer = AllLayers(HighLayerIdx);
        ExpectedHeight = LowLayer.TrueHeight * Para.Zooming;
        ExpectedWidth = LowLayer.TrueWidth * Para.Zooming;
        Scale = OverBigLayer.TrueHeight/ExpectedHeight;
        OriginalScaleSetting = (1.25)^6;
        GauVar = Para.B_GauVar * log(Scale)/log(OriginalScaleSetting);
        %Build Expected High-Res result
        img = OverBigLayer.GridAsHighLayer;
        SubSampleImage = GenerateSubSampleImage(img, ExpectedHeight, ExpectedWidth, Scale, GauVar);

        %fill data
        MappingLayerLow(i).Index = i;
        MappingLayerLow(i).TrueHeight = LowLayer.TrueHeight;
        MappingLayerLow(i).TrueWidth = LowLayer.TrueWidth;
        MappingLayerLow(i).FormatHeight = LowLayer.FormatHeight;
        MappingLayerLow(i).FormatWidth = LowLayer.FormatWidth;
        MappingLayerLow(i).ValidHeight = LowLayer.ValidHeight;
        MappingLayerLow(i).ValidWidth = LowLayer.ValidWidth;
        MappingLayerLow(i).Image = LowLayer.GridAsHighLayer;

        MappingLayerHigh(i).Index = i;
        MappingLayerHigh(i).TrueHeight = ExpectedHeight;
        MappingLayerHigh(i).TrueWidth = ExpectedWidth;
        MappingLayerHigh(i).FormatHeight = ceil(ExpectedHeight);
        MappingLayerHigh(i).FormatWidth = ceil(ExpectedWidth);
        MappingLayerHigh(i).ValidHeight = floor(ExpectedHeight);
        MappingLayerHigh(i).ValidWidth = floor(ExpectedWidth);
        MappingLayerHigh(i).Image = SubSampleImage;
    end
    clear LoadData img SubSampleImage Gauvar OriginalB_GauVar OrignalScaleSetting Scale ExpetectedHeight ExpectedWidth LowLayer OverBigLayer i

    %extract the high-low patch feature
    TempFolder = Para.TempDataFolder;
    SaveName = Para.SaveName;
    LayerPairNum = PairNumber;   %keep the original implementation, in fact, they are repeated.
    LowEdge = Para.LowPatchSize;
    PatchSize = LowEdge^2;
    LowFeatureLength = 4*PatchSize;

    HighEdge = LowEdge * Para.Zooming;
    HighPatchSize = HighEdge^2;
    HighFeatureLength = HighEdge^2;
    ImageDataBase = cell( LayerPairNum,1);
    TotalVectorNum = 0;
    %Compute the TotalVectorNum
    for idx =1:LayerPairNum
        Himg = MappingLayerHigh(idx).Image;
        Limg = MappingLayerLow(idx).Image;
        ImageDataBase{idx}.Himg = Himg;
        ImageDataBase{idx}.Limg = Limg;
        Hrow = MappingLayerHigh(idx).ValidHeight;
        Hcol = MappingLayerHigh(idx).ValidWidth;
        Lrow = MappingLayerLow(idx).ValidHeight;
        Lcol = MappingLayerLow(idx).ValidWidth;
        ImageDataBase{idx}.Hrow = Hrow;
        ImageDataBase{idx}.Hcol = Hcol;
        ImageDataBase{idx}.Lrow = Lrow;
        ImageDataBase{idx}.Lcol = Lcol;
        if (Lrow - 4- (LowEdge-1) < 1) || (Lcol - 4- (LowEdge-1) < 1)
            error('The downsampled image is too small to extract features');
        end
        LPatchNum = (Lrow - 4- (LowEdge-1) )*(Lcol - 4- (LowEdge-1) );    %use intensity feature, the bounday can be fully covered
        TotalVectorNum = TotalVectorNum + LPatchNum;
    end

    %Extract feature: low as gradeint, high as difference from patch mean
    LowFeatureMatrix = zeros( LowFeatureLength + 3, TotalVectorNum );
    HighFeatureMatrix = zeros( HighFeatureLength, TotalVectorNum );
    for i=1:LayerPairNum
        Limg = ImageDataBase{i}.Limg;
        ImageDataBase{i}.LG1 = conv2(Limg , [-1,0,1], 'same');      %here are little problems, for the Matlab-built in conv2. Should I create a own-made conv2?
        ImageDataBase{i}.LG2 = conv2(Limg , [-1,0,1]', 'same');
        ImageDataBase{i}.LG3 = conv2(Limg , [1,0,-2,0,1], 'same');
        ImageDataBase{i}.LG4 = conv2(Limg , [1,0,-2,0,1]', 'same');
    end

    FeatureIdx = 0;
    for ImageIdx = LayerPairNum:-1:1                        %here, L0 has to be put in the first item of ImageDataBase. Otherwise, the cooredinate of the low L0 image will not be in the first part in the LowFeatureMatrix 
        Limg = ImageDataBase{ImageIdx}.Limg;
        Himg = ImageDataBase{ImageIdx}.Himg;
        Lrow = MappingLayerLow(ImageIdx).ValidHeight;
        Lcol = MappingLayerLow(ImageIdx).ValidWidth;
        for i = 3:Lrow-2-(LowEdge-1)
            for j = 3:Lcol-2-(LowEdge-1)
                top = i;
                bottom = i+LowEdge-1;
                left = j;
                right = j+LowEdge-1;
                LG1 = ImageDataBase{ImageIdx}.LG1;
                LG2 = ImageDataBase{ImageIdx}.LG2;
                LG3 = ImageDataBase{ImageIdx}.LG3;
                LG4 = ImageDataBase{ImageIdx}.LG4;
                LPatch1 = LG1(top:bottom,left:right);
                LPatch2 = LG2(top:bottom,left:right);
                LPatch3 = LG3(top:bottom,left:right);
                LPatch4 = LG4(top:bottom,left:right);
                %low-res patches become intensity features
                %Is it reasonable? Hard to say.
                LVector = [reshape(LPatch1, PatchSize , 1);reshape(LPatch2, PatchSize , 1); ...
                    reshape(LPatch3, PatchSize , 1);reshape(LPatch4, PatchSize , 1)];
                FeatureIdx = FeatureIdx + 1;
                LowFeatureMatrix(1:LowFeatureLength,FeatureIdx) = LVector;
                LowFeatureMatrix(LowFeatureLength+1,FeatureIdx) = i;
                LowFeatureMatrix(LowFeatureLength+2,FeatureIdx) = j;
                LowFeatureMatrix(LowFeatureLength+3,FeatureIdx) = ImageIdx;

                hightop = (i-1)*Para.Zooming + 1;
                highleft = (j-1)*Para.Zooming + 1;
                highbottom = hightop + HighEdge -1;
                highright = highleft + HighEdge -1;
                HighPatch = Himg(hightop:highbottom,highleft:highright);
                HighFeature = HighPatch - mean(HighPatch(:));
                HighVector = reshape(HighFeature , HighPatchSize , 1);
                HighFeatureMatrix(1:HighFeatureLength,FeatureIdx) = HighVector;                
            end
        end
    end

    if FeatureIdx ~= TotalVectorNum
        error('error, FeatureIdx should = TotalVectorNum');
    end

    save( fullfile(TempFolder,[SaveName '_FeatureMatrix.mat']),'LowFeatureMatrix','HighFeatureMatrix');

    %Convert To XlXh
    LowVectorLength = LowFeatureLength;
    Xl = LowFeatureMatrix(1:LowVectorLength,:);
    Xh = HighFeatureMatrix;

    dimXl = size(Xl,1);
    dimXh = size(Xh,1);
    CombinedTable = [Xl/sqrt(dimXl);Xh/sqrt(dimXh)];        %for balancing Xl and Xh,

    %mornalize the SortSignal
    SqrtSumOfSquare = sqrt(sum(CombinedTable.^2));
    DataSize = length( SqrtSumOfSquare );
    XlXh = zeros(dimXl+dimXh , DataSize );
    Thd = 1e-5;
    for i=1:DataSize
        if SqrtSumOfSquare(i) < Thd
            XlXh(:,i) = 0;      %Avoid NaN
        else
            XlXh(:,i) = CombinedTable(:,i) ./ SqrtSumOfSquare(i);      %||D_j|| = 1 
        end
    end
    save(fullfile(TempFolder, [SaveName '_XlXh.mat']) , 'XlXh' );   
else
    TempFolder = Para.TempDataFolder;
    SaveName = Para.SaveName;
    LowEdge = Para.LowPatchSize;
    PatchSize = LowEdge^2;
    LowFeatureLength = 4*PatchSize;
    load( fullfile(TempFolder, [SaveName '_FeatureMatrix.mat']),'LowFeatureMatrix','HighFeatureMatrix');
end %end of if Para.bExtractFeature


if Para.bDoCluster
    %Do Cluster
    PureLowFeatureMatrix = LowFeatureMatrix(1:LowFeatureLength,:);      %the last 3 is top, left, and LayerIdx
    SignalNumber = size(LowFeatureMatrix,2);
    KMeansInput = PureLowFeatureMatrix';
    ClusterNum = ceil(SignalNumber / Para.AvgSignalNumPerCluster);
    
    opts = statset('Display' , 'iter');
    disp('Doing kmenas');
    [IDC, C, sumd, D]= kmeans(KMeansInput, ClusterNum , 'emptyaction' , 'drop', 'options' , opts);   %'drop' 'singleton' , 'options',options
    disp('Finish kmeans');
    %do not save this file, it is too big (~2G)
    %if isfield( Para, 'SaveClusterDetail')
    %    if Para.SaveClusterDetail
    %        save( [TempFolder SaveName '_ClusterResult_Avg' num2str(Para.AvgSignalNumPerCluster) '.mat'],'IDC','C','sumd','D', 'ClusterNum');
    %    end
    %end
    
    %confirm the cluster info here (once), instead of repeatedly loading a big file in FindCoef2
    AllMatchArray = false(SignalNumber,ClusterNum);
    for i = 1:SignalNumber
        [B, IX] = sort(D(i,:));
        for j=1:Para.ClusterPerPoint
            AllMatchArray(i,IX(j)) = true;
        end
    end
    save( fullfile(TempFolder, [SaveName '_ClusterResultAllMatchArray_Avg' num2str(Para.AvgSignalNumPerCluster) '.mat']), 'AllMatchArray', 'ClusterNum');    
    clear IDC C sumd D ClusterNum AllMatchArray ClusterNum B IX opts PureLowFeatureMatrix SignalNumber KMeansInput
end
    
for Iter=Para.IterNumStart:Para.IterNumEnd
    if Iter == 1
        CreateInitialDictionary( Para );
    end
    CoefMatrix = FindCoef3( Para ,Iter );
    ReconstructHigh3(Para , img_input_y, IQLayer, Iter);        %ReconsructHigh3 = ReconsructHigh + back-projection
    UpdateDictionary( Para ,Iter);
end