FER: Facial Expression Recognition

الگوریتم سیفت: ( scale invariant feature transform )

در مبحث به دست آوردن نقاط ویژه که کاربرد فراوانی در انطباق تصویر دارد ، به دست آوردن نقاط ویژه هریس (harris point  )

از پایه ای ترین روش ها می باشد با این وجود عدم کارایی این روش و روش های مشابه در مقیاس های مختلف ( scale ) از تصویر استفاده از این روش ها رو به چالش می کشد.

با هدف رفع این مشکل ویژگی های SIFT ( Scale Invariant Feature Transfirm)  معرفی شدند. در این روش تصویر به جای محاسبه و استفاده از لاپلاسین گوسین LoG  از تفاضل گوسی DoG استفاده می شود و به جای استفاده از سیگماهای مختلف در فیلتر گوسی ، تفاضل گوسی روی اکتاو octave  مختلف از هرم تصویر اعمال می شود تا نقاط ویژه به دست آید .

در این بین برای حذف نقاط نا مناسب از یک سطح آستانه استفاده می شود که اگر حاصل تفاضل از این مقدار کمتر بود یعنی نقطه قدرتمندی نیست و در نظر گرفته نمی شود. از سوی دیگر لبه ها همچنین برای تشخیص نقاط مناسب مشکل ایجاد می کنند برای منظور در ماسک 2*2 تعریف شده مقادیر ویژه محاسبه می شود اگر دو تفاوت دو مقدار ویژه بیش از یک آستانه مشخص بود این نقطه لبه بوده و در نظر گرفته نمی شود.

همچنین برای هر نقطه زاویه یا جهت تعیین می شود که معمولا 36 جهت در نظر گرفته می شود که بر اساس گرادیان محاسبه می شود.در نهایت برای هر نقطه همسایگی 16*16 به صورت 16 بلوک 4*4 بررسی می شود که در هر بلوک هیستوگرام جهت ها به صورت 8 ستونی محاسبه می شوند پس در انتها به ازائ هر نقطه 128 مقدار بردار توصیف گر نقطه را تشکیل می دهند. در تصویر زیر به ترتیت از چپ به راست تصویر ورودی ، تصویر وسط تصویرنقاط استخراج شده و تصویر سمت راست جهت غالب در هر نقطه را نشان می دهد


نمونه ای از اجرای  الگوریتم sift   به صورت زیر است : 

%% Initialization

% Here, the x-axis correspond to coloum, while the y-axis correspond to row

clear all; close all;

img = imread('cameraman.tif');

[~,~,ColorChannel] = size(img);

if ColorChannel > 1

    img = rgb2gray(img);

end

img = double(img)/255;

ScaleSpaceNum = 3; % number of scale space intervals

SigmaOrigin = 2^0.5; % default sigma

ScaleFactor = 2^(1/ScaleSpaceNum);

StackNum = ScaleSpaceNum + 3; % number of stacks = number of scale space intervals + 3

OctaveNum = 3;

GaussianFilterSize = 21;

OctaveImage = {OctaveNum,StackNum}; % save the Gaussian-filtered results of image

OctaveImageDiff = {OctaveNum,StackNum-1}; % save the Difference of Gaussian-filtered results of image



%% Gaussian Convolution of Images in Each Octave

ImgOctave = cell(OctaveNum,1);

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1); % when up to a new octave, double the sigma

    ImgOctave{Octave} = imresize(img, 2^(-(Octave-1)));

    for s = 1:StackNum

        SigmaScale = Sigma * ScaleFactor^(s-2);

        % calculate Guassian kernel

        GaussianFilter = fspecial('gaussian',[GaussianFilterSize,GaussianFilterSize],SigmaScale);

        % do convolution with Gassian kernel

        OctaveImage{Octave,s} = imfilter(ImgOctave{Octave}, GaussianFilter,'symmetric');

    end  

end

% calculate difference of Gaussian (in original paper eq.1)

for Octave = 1:OctaveNum

    for s = 1:StackNum-1

        OctaveImageDiff{Octave,s} = OctaveImage{Octave,s+1} - OctaveImage{Octave,s};

    end

end

%% Find the local minima and maxima between stacks

DiffMinMaxMap = cell(OctaveNum,StackNum-3);

for Octave = 1:OctaveNum

    for s = 2:size(OctaveImageDiff,2)-1

        CompareDiffImg = zeros(size(OctaveImage{Octave,1},1)-2,size(OctaveImage{Octave,1},2)-2,27);

        indx = 0; % 3rd dimension indx for CompareDiffImg

        for s2 = s-1:s+1

            for k = 1:9

                [i,j] = ind2sub([3,3],k);

                CompareDiffImg(:,:,indx+k) = OctaveImageDiff{Octave,s2}(i:end-3+i,j:end-3+j);

            end

            indx = indx + 9;

        end

        [~,MinMap] = min(CompareDiffImg,[],3);

        MinMap = (MinMap == 14);

        [~,MaxMap] = max(CompareDiffImg,[],3);

        MaxMap = (MaxMap == 14);

        DiffMinMaxMap{Octave,s-1} = MinMap + MaxMap; % the center indx is 9 + 5 = 14

    end

end


%% Accurate keypoint localization

UnstableExtremaThreshold = 0.03; % set the threshold as 0.03 (the same as the original paper)

DiffMinMaxMapAccurate = DiffMinMaxMap;

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1);

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2)

        Map = DiffMinMaxMap{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            if (Row <= 1) || (Row >= size(Map,1)) || (Col <= 1) || (Col >= size(Map,2))

                DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard out of matrix boundary

                continue;

            end

            RowDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col);

            ColDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col+1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col-1);

            ScaleDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+2}(Row,Col) - OctaveImageDiff{Octave,DiffMinMaxMapIndx}(Row,Col);

            offset = [2; 2; Sigma * ScaleFactor^(DiffMinMaxMapIndx) - Sigma * ScaleFactor^(DiffMinMaxMapIndx-2)];

            DxHat = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + 0.5 * ([RowDiff,ColDiff,ScaleDiff] * offset);

            if abs(DxHat) < UnstableExtremaThreshold

                DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard unstable extrema

            end

        end

    end

end

%% Eliminating edge responses

gamma = 10; % set the threshold gamma as 10 (the same as the original paper)

DetermineThreshold = (gamma + 1)^2 / gamma;

DiffMinMaxMapNoEdge = DiffMinMaxMapAccurate;

for Octave = 1:OctaveNum

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2) 

        Map = DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            if (Row <= 1) || (Row >= size(Map,1)) || (Col <= 1) || (Col >= size(Map,2))

                DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard out of matrix boundary

                continue;

            end

            Dyy = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col) - 2*OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col);

            Dxx = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col+1) - 2*OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col-1);

            Dxy = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col+1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col-1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col+1) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col-1);

            TrH = Dxx + Dyy;

            DetH = Dxx*Dyy - Dxy^2;

            if ((TrH^2 / DetH) >= DetermineThreshold)

                DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard unstable extrema

            end

        end

    end

end


%% SIFT feature descriptors generation

% the patch size for dominant orientation calculation is 9;

% the patch size for feature transformation is 16;

DominantOrientation = cell(size(DiffMinMaxMap));

SIFT = cell(size(DiffMinMaxMap));

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1); % when up to a new octave, double the sigma

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2)

        stack = DiffMinMaxMapIndx+1;

        SigmaScale = Sigma * ScaleFactor^(stack-2);

        GaussianSmoothedImage = OctaveImage{Octave,stack};

        Map = DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        DomOri = zeros(length(SSCindx),2,36); % first column is for sita, second column is for magnitude

        sift = zeros(length(SSCindx),128,36);

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            Row = Row+1; Col = Col+1; % offset = [1,1];

            if (Row <= 10) || (Row >= size(GaussianSmoothedImage,1)-8) || (Col <= 10) || (Col >= size(GaussianSmoothedImage,2)-8)

                % skip if out of boundary

                continue;

            end

            % magnitude and sita of sample points in the patch

            mag = ((GaussianSmoothedImage(Row-8:Row+7,Col-7:Col+8) - GaussianSmoothedImage(Row-8:Row+7,Col-9:Col+6)).^2 + (GaussianSmoothedImage(Row-7:Row+8,Col-8:Col+7) - GaussianSmoothedImage(Row-9:Row+6,Col-8:Col+7)).^2).^0.5;

            sita = atan2((GaussianSmoothedImage(Row-7:Row+8,Col-8:Col+7) - GaussianSmoothedImage(Row-9:Row+6,Col-8:Col+7)),(GaussianSmoothedImage(Row-8:Row+7,Col-7:Col+8) - GaussianSmoothedImage(Row-8:Row+7,Col-9:Col+6)));

            sita = mod(sita + 2*pi, 2*pi); % the range of atan2 function is -pi~pi, map it to 0~2*pi

            % Dominant orientation calculation

            GaussianFilter = fspecial('gaussian',[9,9],1.5*SigmaScale);

            Dmag = mag(5:13,5:13).* GaussianFilter; % magnitude is weighted by gaussian filter

            sitaquantize = mod(sita(5:13,5:13) + pi/36,2*pi);

            sitaquantize = floor(sitaquantize / (2*pi/36));

            sitabin = zeros(36,1);

            for patchindx = 1:9^2

                sitabin(sitaquantize(patchindx)+1) = sitabin(sitaquantize(patchindx)+1) + Dmag(patchindx);

            end

            maxsitabin = max(sitabin);

            DominantOriBin = find((sitabin / maxsitabin) >= 0.8); % duplicate the feature when non-maximum magnitude of orientation is also big

            DominantOriSita = (DominantOriBin-1)*(2*pi/36);

            DomOri(ssc,1,1:length(DominantOriSita)) = DominantOriSita;

            DomOri(ssc,2,1:length(DominantOriSita)) = sitabin((sitabin / maxsitabin) >= 0.8);

            % SIFT feature calculation

            for DomOriIndx = 1:length(DominantOriSita)

                GaussianFilter = fspecial('gaussian',[16,16],8);

                Smag = mag.* GaussianFilter;

                sitaquantize = mod(sita - DominantOriSita(DomOriIndx) + pi/8 + 2*pi,2*pi);

                sitaquantize = floor(sitaquantize / (2*pi/8));

                sitabin = zeros(8,4,4);

                for patchindx = 1:16^2

                    [row,col] = ind2sub([16,16], patchindx);

                    row = floor((row-1)/4)+1;

                    col = floor((col-1)/4)+1;

                    sitabin(sitaquantize(patchindx)+1,row,col) = sitabin(sitaquantize(patchindx)+1,row,col) + Smag(patchindx);

                end

                sitabin = bsxfun(@times,sitabin,sum(sitabin.^2,1).^(-0.5)); % normalize the vector

                sitabin(sitabin > 0.2) = 0.2; % threshold the maximum value as 0.2

                sitabin = bsxfun(@times,sitabin,sum(sitabin.^2,1).^(-0.5)); % renormalize

                sift(ssc,:,DomOriIndx) = sitabin(:);

            end

        end

        DominantOrientation{Octave,DiffMinMaxMapIndx} = DomOri;

        SIFT{Octave,DiffMinMaxMapIndx} = sift;

    end 

end


%% Show results (this part of code is just for verification, and is a mess...)

% close all;

% imgtemp = img*0.5;

% imgtemp(2:end-1,2:end-1) = (imgtemp(2:end-1,2:end-1) + (DiffMinMaxMap{1,1}+DiffMinMaxMap{1,2}+DiffMinMaxMap{1,3}))*255;

% figure,imshow(uint8(imgtemp));

% imgtemp2 = img*0.5;

% imgtemp2(2:end-1,2:end-1) = (imgtemp2(2:end-1,2:end-1) + (DiffMinMaxMapAccurate{1,1}+DiffMinMaxMapAccurate{1,2}+DiffMinMaxMapAccurate{1,3}))*255;

% figure,imshow(uint8(imgtemp2));

% imgtemp3 = img*0.5;

% imgtemp3(2:end-1,2:end-1) = (imgtemp3(2:end-1,2:end-1) + (DiffMinMaxMapNoEdge{1,1}+DiffMinMaxMapNoEdge{1,2}+DiffMinMaxMapNoEdge{1,3}))*255;

% figure,imshow(uint8(imgtemp3));

imgtemp4 = img*0.5;

figure,imshow(uint8(imgtemp4*255));

for i = 1:3

for j = 1:3

for k = 1:2

if isempty(DominantOrientation{i,j})

continue

end

DomOri = DominantOrientation{i,j}(:,:,k);

DOx = cos(DomOri(:,1)).*DomOri(:,2)*200;

DOy = sin(DomOri(:,1)).*DomOri(:,2)*200;

Map = DiffMinMaxMapNoEdge{i,j};

SSCindx = find(Map); % Scale-Space-Corner Index

for ssc = 1:length(SSCindx)

    hold on;

    [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

    Row = Row*2^(i-1)+1; Col = Col*2^(i-1)+1; % offset = [1,1];

    if (Row > size(img,1)) || (Col > size(img,2))

        continue

    end

    quiver(Col,Row,DOx(ssc),DOy(ssc),'r')

end

end

end

end

hold off;

%

imgtemp5 = img*0.5;

for i = 1:3

for j = 1:3

for k = 1:2

if isempty(DominantOrientation{i,j})

continue

end

Map = DiffMinMaxMapNoEdge{i,j};

SSCindx = find(Map); % Scale-Space-Corner Index

for ssc = 1:length(SSCindx)

    [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

    Row = Row*2^(i-1)+1; Col = Col*2^(i-1)+1; % offset = [1,1];

    if (Row > size(img,1)) || (Col > size(img,2))

        continue

    end

    imgtemp5(Row,Col) = 1;

end

end

end

end

figure,imshow(uint8(imgtemp5*255));

K نزدیک ترین همسایه (K-nearest neighbour)

الگوریتم نزدیکترین همسایگی برای طبقه بندی یا classification  رگرسیون regression  کاربرد دارد . این راهکار در فاز آموزش شامل تعداد نقاط آموزشی با کلاس یا برچسب مشخص می باشد و پس از آن برای نمونه های جدید در فاز تست ، کلاس خروجی بر اساس آرای k  تا از نزدیکترین همسایه هایش مشخص می شود که اگر k=1  باشد کلاس نمونه بر اساس کلاس نزدیکترین همسایه تعیین می شود.

در متلب برای ایجاد یک مدل طبقه بندی knn  تابع مناسب تعریف شده و می توان از آن استفاده نمود به طور مثال قطعه کد زیر :

  load fisheriris

X = meas;

Y = species;

Mdl = fitcknn(X,Y,'NumNeighbors',5)

در این مثال از دیتاست تعریف شده fisheriris استفاده شده است که شامل داده های گل ها در متغیر X,meas ورودی ها و Y,species  برچسب ها یا انواع گل ها می باشد برای ساخت مدل طبقه بندی از تابع fitcknn  استفاده شده است که تعداد همسایه ها نیز 5 در نظر گرفته شده است .

قطعه کد زیر نحوه استفاده و خروجی را برای مدل ساخته شده با استفاده از تابع predict  نشان می دهد.

predict(Mdl,X(2,:))

 

ans =

 

  cell

 

    'setosa'

چنانچه ملاحظه می شود مدل به درستی کلاس نمونه دوم را که ‘setosa’ می باشد را تشخیص داده است.

طبقه بندی مبتنی بر نمایش پراکنده (Sparse Representation-based Classification)

SRC بر مبنای حس کردن فشرده (CS) توسعه یافته است. اصل روش SRC بر مبنای این فرض است که کل مجموعه نمونه های آموزشی برای تشکیل یک دیکشنری استفاده می شوند و بنابراین، مسئله ی طبقه بندی به صورت مسئله ی جست و جوی متمایز کننده ی نمایش پراکنده ی نمونه تست به صورت ترکیب خطی نمونه های آموزشی با حل مسئله بهینه سازی معیار L1 در نظر گرفته می شود. در حالت رسمی، برای نمونه های آموزشی یک کلاس، این فرض می تواند به صورت زیر فرموله شود:

که در آن  نشان دهنده ی نمونه تست K امین کلاس است، نشان دهنده ی i امین نمونه ی آموزشی K امین کلاس است،نشان دهنده ی وزن متناظر وزن است و نشان دهنده ی خطای تقریب است.

برای نمونه های آموزشی از تمام C کلاس شی، معادله بالا می تواند به صورت زیر بازنویسی شود:

در فرم ماتریس، برابر است با:

با شرط:

برای دستیابی به بردار وزن ،مسئله ی کمینه سازی معیار L1 زیر باید حل شود:

که یک مسئله ی بهینه سازی محدب است و می تواند با روش برنامه نویسی درجه دوم حل شود. زمانی که راه حل پراکنده یارائه می شود، روند طبقه بندی SRC به صورت زیر خلاصه می شود:

گام 1: حل مسئله ی معادله کمینه سازی معیار L1.

گام 2: برای هر کلاس i، محاسبه ی باقیمانده های بین نمونه ی بازسازی شده ی و  تست مورد نظر

گام 3: برچسب کلاس نمونه تست مشخص با استفاده از قاعده زیر تعیین می شود: شناسایی 

در کارهای قبلی، کارایی SRC هنگام طبقه بندی تصاویر حالت چهره واضح یا مات مورد بررسی قرار گرفته است و پی برده شده است که SRC، کارایی بهتر و قدرت بالاتری در مقایسه با نزدیک ترین همسایه (NN)، نزدیک ترین زیر فضا(NS) و SVM داشته است. گروه تحقیقاتی محمدی، یک دیکشنری مبتنی بر PCA با ساخت نمایش پراکنده  و طبقه بندی حالت های چهره جهانی ارائه کرده است. در جزییات، تصاویر چهره های گویای هر فرد ابتدا از تصویر چهره ی خنثی همان فرد کسر شدند. سپس، روش PCA برای این تصاویر متفاوت جهت مدل سازی واریانس های درون هر کلاس حالت چهره به کار گرفته شد. مؤلفه های اصلی یاد گرفته شده به عنوان اجزای سازنده ی دیکشنری به کار گرفته شدند. در نهایت، برای طبقه بندی، تصویر تست مورد نظر به طور پراکنده به صورت ترکیب خطی مؤلفه های اصلی شش حالت چهره پایه نمایش داده شد. گروه تحقیقاتی اویانگ اخیرا FER دقیق و قوی را با ترکیب طبقه بندی کننده های مبتنی بر نمایش پراکنده چند گانه توسعه داده است، یعنی ترکیب HOGCSRC و LBPCSRC.

<< 1 2 3 4 5 ... 8 >>