الگوریتمSURF(Speed Up Robust Feature)

استفاده از ویژگی های SIFT زمانبر بوده و پیچیدگی  محاسباتی بالایی را در بر داشت از این رو محققان روش دیگری به اسم SURF (Speed Up Robust Feature) ارائه نمودند. در این روش به جای استفاده از DOG یا تفاضل گوسی از باکسی استفاده می شود که با تصویر کاتوایلو می شود و برای محاسبه ی حاصل کانولوشن از تصویر انتگرالی که یک بار در ابتدا محاسبه شده است استفاده می شود. نکته حائز اهمیت در این روش امکان محاسبه حاصل کانولوشن باکس ها با ابعاد مختلف به صورت موازی می باشد.

در این روش همچنین برای محاسبه جهت، از ویولت استفاده می شود. استفاده از ویولت باعث می شود تا نقاطی که نیاز نیست جهت محاسبه صورت نگیرد که این امر تاثیر شایانی در افزایش صورت الگوریتم دارد.

به طور مثال در تصویر زیر سمت چپ تصویر ورودی و در سمت راست 10 نقطه توسط SURF به دست آمده مشخص شده است.

I = imread('cameraman.tif');
points = detectSURFFeatures(I);
%Display locations of interest in image.

imshow(I); hold on;
plot(points.selectStrongest(10));

ماشین بردار پشتیبان(SVM)

استفاده از ماشین های بردار پشتیبان بسیار در زمینه ی دسته بندی بسیار پر کاربرد است و دقت خوبی در طبقه بندی فراهم می نماید. در این روش در حالت دو بعدی یک خط و در حالت کلی یک ابر صفحه بدست آورده می شود که دو کلاس را از یکدیگر جدا می نماید و این ابر صفحه یا خط به گونه ای انتخاب می شود که بیشترین حاشیه یا margin را نسبت به داده های موجود در هر دو کلاس داشته باشد.

برای استفاده از SVM به صورت چند کلاسه از تعدادی SVM دو کلاس استفاده می شود. به طور مثال در تابع fitcecoc به ازای k کلاس خروجی، از k(k-1)/2  ماشین بردار پشتیبان استفاده می شود. استفاده از این تابع کار با SVM را در حالت چند کلاسه آسان می نماید. در ادامه مثالی از استفاده از SVM در حالت چند کلاسه بیان می شود.

 load fisheriris

X = meas

Y = species

 Mdl = fitcecoc(X,Y)

Mdl =

  ClassificationECOC

             ResponseName: 'Y'

    CategoricalPredictors: []

               ClassNames: {'setosa'  'versicolor'  'virginica'}

           ScoreTransform: 'none'

           BinaryLearners: {3×1 cell}

               CodingName: 'onevsone'

  Properties, Methods

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

>> output=predict(Mdl,X(90,:))

output =

  cell

    'versicolor'

شبکه عصبی مصنوعی (ANN)

شبکه های عصبی مصنوعی ( ANN) ابزار بسیار قوی در زمینه ی هوش مصنوعی است که در طبقه بندی، درون یابی (پیش بینی) کاربرد فراوانی دارد.

بنا بر کاربردهای مختلف، انواع مختلفی از شبکه های عصبی معرفی شده است که شبکه عصبی پرسپترون چند لایه (mlp) شبکه بسیار قوی از دسته شبکه های عصبی پیش خور (feed forward) است که سال هاست در زمینه طبقه بندی مورد استفاده قرار می گیرد. این شبکه عصبی معمولا با استفاده از یک لایه یا دو لایه از نرون های مخفی با توابع فعالیت غیر خطی، قادر به طبقه بندی مسائل غیر خطی با دقت بالایی است. تعداد نرون ها در لایه های مخفی بسته به پیچیدگی مسئله و ورودی های شبکه های عصبی دارد و معمولا از یک یا دو لایه مخفی استفاده می شود و استفاده از لایه های بیشتر در مبحث یادگیری عمیق و معماری های مرتب جای می گیرند.

تابع newff از نسخه های بسیار قبل تر در نرم افزار متلب تعریف شده و یک شبکه mlp ایجاد می کند. برای استفاده از این تابع نیاز است ورودی های، هدف ها و تعدادلایه های مخفی و تابع فعالیت نرون ها و همچنین روش آموزشی بیان گردد. در صورتی که توابع فعالیت و روش آموزش بیان نشود، متلب به صورت پیش فرض از الگوریتم لونبرگ - مارکوارت و توابع فعالیت سیگموید دو قطبی استفاده می کند در ادامه نمونه ای از استفاده این شبکه ها در متلب آورده شده است.

>> [input,target] = simplefit_dataset;

>> net=newff(input,target,10,{'tansig','purelin'},'trainlm');

>> net=train(net,input,target);

پس از اجرای دستورات بالا پنجره آموزش شبکه عصبی باز شده و مولفه های آموزش چنانچه در شکل زیر مشخص است مشاهده می شود که برای مثال بالا شبکه به دقت 4.34e-4 یا حدود چهار ده هزارم در مقیاس mse یا میانگین مربعات خطا دست یافته است.

پس از آموزش شبکه عصبی برای استفاده می توان شبکه و ورودی را به تابع sim ارائه کرد و خروجی شبکه عصبی را به دست آورد. چنانچه مشاهده می شود به طور مثال برای 14 نمونه شبکه دقیقا مقدار هدف ایده آل را بدست می آورد.

((output=sim(net,input(14

=output 

    9.1432

( target(14

=ans 

    9.1432

الگوریتم سیفت: ( 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));