Sauliaus Sinkevičiaus blogas

Klasifikavimo uždaviniui spręsti yra daugybė klasifikatorių, o mes palengva šį bei tą išmėginsim, pažiūrėsim kaip jie veikia ir kaip atrodo iš arčiau.

Šiandien mūsų taikinyje yra kvadratinis Gauso klasifikatorius. Teoriją apie klasifikatorių įkėliau atskirai, bet ji angliškai ir gana skurdi. Šis klasifikatorius leidžia rasti ribą tarp dvejų klasių, tačiau norint gero tikslumo, reikia pasiekti, kad požymiai turėtų normalinį pasiskirstymą. Naudojami duomenys yra čia, tik du požymiai ir dvi klasės. Programos kodas:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

clc;
close all;
 
% Užkraunami duomenys
load Apple.txt
data = Apple;
clear Apple;
 
% Braižomas grafikas
y = data(:,3) == 0;
d1 = data(y,1:2)';
d2 = data(~y,1:2)';
figure(1); clf; hold on
plot(d1(1,:), d1(2,:), 'r+');
plot(d2(1,:), d2(2,:), 'g+');
axis('equal');
 
% Parametrai
fromCol = 1; % nuo kurio stulpelio prasideda požymiai, imtinai
features = 2; % požymių kiekis
till = fromCol + features - 1; % paskutinis požymių stulpelis
classCol = 3; % stulpelis, kuriame yra klasės numeris
 
% Tikslumo skaičiavimas su "10-fold cross validation"
tic;
accuracy = mycrossvalidation('Gaussian', data, classCol, fromCol, features);
fprintf('Tikslumas = %3.2f, laikas = %3.2f s \n', accuracy, toc);
 
% Piešiamas klasifikatoriaus diskriminantinės funkcijos kontūras
cols = fromCol + features - 1; % stulpelių kiekis duomenyse
k = data(:, classCol) == 0;
c1 = data(k, fromCol : cols)'; % apmokymo duomenų pirma klasė
c2 = data(~k, fromCol : cols)'; % apmokymo duomenų antra klasė
d = -2 : .05 : 2;
[X, Y] = meshgrid(d, d);
xx = 1;
yy = 1;
Z = zeros(size(X));
for x = d
    for y = d
        Z(yy, xx) = quadratic_gaussian_classifier([x y]', c1, c2);
        yy = yy + 1;
    end;
    xx = xx + 1;
    yy = 1;
end;
contour(X, Y, Z, [0,0], 'k');

Naudojama „10-fold cross validation“ validacija, kuri leidžia patikrinti klasifikatoriaus efektyvumą su skirtingomis duomenų imtimis:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

function [accuracy] = mycrossvalidation(type, data, classCol, fromCol, features)
el = size(data, 1); % kiek duomenų yra iš viso
p = 90; % kiek duomenų skirti apmokymui procentais
q = int16(size(data, 1) * p / 100); % kiek duomenų skirti apmokymui
 
total = 0; % bus skaičiuojami teisingi atpažinimai
step = el - q;
for w = 1 : step : el % 10-fold cross validation
    % paskaičiuojami nuo (from) iki (till) intervalas testavimui
    from = w;
    till = w + el - q - 1;
    if till > el
        till = el;
    end;
    % atrenkama apmokymo ir testavimo dalis
    testing = data(from : till, :);
    training = data;
    training(from : till, :) = [];
    % tikrinamas konkretus klasifikatorius, šiuo atveju turime tik vieną
    if strcmp(type, 'Gaussian') == 1
        [count, accuracy] = mygaussian(training, testing, classCol, fromCol, features);
    end;
    total = total + count; % sumuojamas teisingai atpažintas duomenų kiekis
end;  
 
accuracy = total * 100 / el;

Skaičiuojamos tikimybės, kuriai klasei priklauso konkretūs atvejai iš testavimo duomenų imties:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

function [count, accuracy] = mygaussian(training, testing, classCol, fromCol, features)
cols = fromCol + features - 1; % stulpelių kiekis duomenyse
 
k = training(:, classCol) == 0;
c1 = training(k, fromCol : cols)'; % apmokymo duomenų pirma klasė
c2 = training(~k, fromCol : cols)'; % apmokymo duomenų antra klasė
 
count = 0;
el = size(testing, 1);
for z = 1 : el % tikrinam su visais testavimo duomenimis
    x = testing(z, fromCol : cols)';
    class = quadratic_gaussian_classifier(x, c1, c2);
    if class == testing(z, classCol) % ar apskaičiuota klasė atitinka tikrą
        count = count + 1;
    end;
end;
 
accuracy = count * 100 / el;

Ir galiausiai kvadratinio Gauso klasifikatoriaus realizacija pagal 2.20 formulę:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

function [class] = quadratic_gaussian_classifier(x, c1, c2)
covM1 = cov(c1');
covM2 = cov(c2');
invM1 = inv(covM1);
invM2 = inv(covM2);
piko1 = mean(c1, 2);
piko2 = mean(c2, 2);
pc1 = size(c1, 2);
pc2 = size(c2, 2);
A = 0.5 * (invM1 - invM2); % 2.21 formulė
beta = invM2 * piko2 - invM1 * piko1; % 2.22 formulė
gama = (piko1' * invM1 * piko1 - piko2' * invM2 * piko2) * 0.5 ...
    + log(det(covM1) / det(covM2)) - log(pc1 / pc2); % 2.23 formulė
res = x' * A * x + beta' * x + gama; % 2.20 formulė
if res < 0
    class = 0;
elseif res > 0
    class = 1;
else
    class = -1; % neapsisprendęs
end;

Skaičiavimų rezultatas:

1

Tikslumas = 97.68, laikas = 1.14 s

Diskriminantinės funkcijos kontūras parodytas paveikslėlyje, taškų braižymui paimta nedaug, todėl kreivė laužyta, tačiau tai yra parabolė:

97,68 % tikslumas yra gana aukštas. Šį klasifikatorių bandžiau su duomenimis, kurie turi 30 požymių ir gavau 95,96 % tikslumą, o atmetus mažiau reikšmingus požymius ir palikus jų tik 8 bei kai kuriuos iš jų normalizavus, pavyko pasiekti 97,01 % tikslumą. Vadinasi klasifikatorius yra gana tikslus, greitas ir patogus. Tačiau jis tinka tik gana paprastiems uždaviniams.

Tyrinėjimai klasifikatorius, klasifikavimas, quadratic Gaussian qlassifier

Nuotraukoj yra mano pagrindinis, esminis ir svarbiausias įrankis mokytis anglų kalbos. Rimtai. Šio daikto efektyvumas stulbinantis. Pačias ausines gavau kaip priedą prie telefono, tačiau jos buvo baisiai nepatogios, krisdavo iš ausų, o laidai skirtingų ilgių, todėl patobulinau. Turėtų matytis kaip viela susukti laidai, kad labai nesimakaluotų, o ausinių galai neilgi, kad po smakru abu laidai susieitų. Neveržia, nekrenta, patogu.

O įrašas apie visai ką kitką. Tai bus mano poziciją dėl A.J. Hoge pasisakymo:

A.J. Hoge stilius man nepatinka, tačiau žaviuosi tuo ką jis daro. Jis inovatoriškas, ieško naujų kelių ir juos randa. Matyt jis vienas iš pirmųjų, kuris bado pirštu, kad švietimo sistema kalbų mokyme (tiesa sakant ir ne tik), nuėjo pernelyg neefektyviu keliu. Žinau daugybę puikių žmonių, labai gabių, kurie mokyklos ar universiteto iškankinti, nepasitiki savim, galvoja, kad yra kalboms negabūs, kad kažkas negerai su jais. Tačiau tiesa gerokai paprastesnė, jie paprasčiausiai suklaidinti ir palikti ant ledo… Jie ir toliau griebiasi to vienintelio kelio, to tikslo, kuris vadinasi „išmokti“. Hoge teisus, tobulai kalbėti neįmanoma. Pridėčiau, kad ir išmokti užsienio kalbos neįmanoma. Jei po teisybei, jūs nemokat ir lietuvių kalbos. Niekas nemoka lietuvių kalbos. Ir niekada nemokės. Pabandykit atsiversti žodyną. Kaip greitai rasite dar nežinomų žodžių? Pabandykit parašyt rašinį ir nunešti kalbininkui, ne tik rašinys raudonuos nuo pabraukimų ties klaidomis, bet ir jūsų veidas. Pabandykit nuvažiuoti į atokesnes Lietuvos vietas, pakalbinti vietinius, įdomu kaip seksis susikalbėti. Pabandykit paskaityti įstatymus, jei pavyzdžiui esat meno žmogus, arba paskaityti apie madas jei esat juristas. Greitai papulsit ne tik į nesuvokiamų sąvokų pinkles, bet ir apskritai į kitą pasaulį. Pabandykit paskaityt laikraščius, kurie yra kelių dešimčių metų senumo, pamatysit, kad lietuvių kalba buvo kitokia.

Tai gal vis gi lietuvių kalbos jūs nemokat? Man atrodo, kad nemokat. Esu tikras, kad visų lietuvių kalbos žodžių nežinot, visų kalbos taisyklių nežinot, atsirastų daugybę lietuvių, kurių nesuprastumėt, ir galybę sričių, kurių įprastos, kasdieninės, kalbos nesuprastumėt. Todėl nėra prasmės net ir statytis sau tikslo „išmokti“. „Išmokti“ kalbą yra utopinis tikslas.

Pabandykim pažiūrėti iš kitos pusės. Bet juk lietuvių kalbą naudojat? Teisingai? Kas įdomu – perskaitot, pažiūrit įdomius filmus, galbūt serialus. Pasikalbat su žmonėmis, kai to reikia ar norit. Kitaip tariant lietuvių kalbos nemokat, bet ją naudojat. Todėl, kai kyla noras vėl griebtis utopinio tikslo – „išmokti“ anglų, prancūzų, vokiečių ar kinų, nesvarbu ką, galima bandyti tą tikslą pervadinti į „naudoti“. Ir tada savęs klausti, o kaip aš galėčiau tą kalbą naudoti?

• gal yra serialų ta kalba, kuriuos norėčiau pamatyti?
• gal yra knygų kurias norėčiau perskaityti originalo kalba?
• gal turiu hobi, apie kurį yra puikios vaizdo-garso-teksto pavidalu užsienio kalba?
• gal kažkas youtube‘ina apie mano dominančias sritis ar hobi?
• gal man patinka žinoti kas vyksta kitose šalyse ir domintų pasiklausyti jų žiniasklaidos?
• gal norėtum tobulėti kaip kulinaras, bet geriausios laidos užsienio kalba?
• gal visada norėjai išmokti piešti ir gali susirasti pamokų ta kalba?
• gal yra kokie nors mmorpg žaidimas, kuriame norit išbandyt komandinį darbą?
• gal yra forumai apie tavo karjerą, hobi ar šiaip gyvenimą?
• ir t.t.

Kai pradedi domėtis, atrodo neįtikėtina, kiek daug Žmogus jau yra sukūręs ir dar mažiau neįtikėtina, kad neįmanoma susirasti kas kiekvienam individualiai labiausiai tiktų ir patiktų.

Užtenka pakeisti tikslą iš „kaip išmokti?“ į „kaip naudoti?“ ir prasideda svarbiausias etapas kalbos tobulėjime, kuriame pats žmogus pradeda ieškoti to kas jam labiausiai tiktų. Jis nebesitiki, kad tai padarys kas nors už jį, nebesirašo į kursus, nebefantazuoja pagyventi svetimoje šalyje, nes neva taip geriausiai išmokstama. Pagaliau jis elgiasi savarankiškai. Pats atsirenka sritį, sudėtingumą, apimtį, pats įvertina, kad jam tikrai bus įdomu ir, kad jis asmeniškai būtinai gaus kažką daugiau nei užsienio kalbos įgūdžius. Galbūt geras emocijas, galbūt žinių, galbūt puikių pražystamų užsienyje ar įgūdžių visai kitose srityse, nes jis pats susiras pačią geriausią ir tinkamiausią medžiagą sau.

Šiuo metu naudoju tris kalbas: lietuvių, anglų, rusų. Aš jų nemoku ir nesimokau, bet naudoju, o kai naudoju, neišvengiamai tobulėju.

Tyrinėjimai išmokti kalbą, kaip mokytis, kalbų mokymąsis, mokytis kalbos, užsienio kalbos

KNN (k nearest neighbors) metodas klasifikuoja tikrindamas kokioms klasėms priklauso artimiausi K kaimynai ir priskiria tai klasei, kuriai priklauso daugelis K kaimynų. Įėjimo vektorius gali turėti bet kokį požymių kiekį, pvz.: 1, 2, 10 ir daugiau. Euklido atstumą naudosime skaičiuojant atstumą tarp kaimynų, kurį lengva paskaičiuoti net ir daug dimensijų turinčiai erdvei. Išbandysime pasirašyti savo KNN funkciją bei palyginti rezultatą su Matlab knnclassify funkcija naudojant dvejų dimensijų duomenų aibę. Abejų klasių duomenys (0 – raudona, 1 – žalia) parodyti žemiau:

Eksperimentas labai paprastas. Užkrausime duomenis, kurie turi du požymius, o pačių duomenų yra kiek daugiau nei tūkstantis. Apmokymui imsime dešimt procentų duomenų atsitiktine tvarka, o tikslumą tikrinsime su likusiais.

Programa:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

clc;
close all;
 
load apple.txt
data = apple;
% duomenys išmaišomi atsitiktine tvarka, nes kas kartą norime
% naudoti vis kitą imti apmokymui, o likutį testavimui
data = data(randperm(size(data, 1)), :);
clear apple;
 
% atvaizduojam abi klases
figure(1); clf; hold on
l = data(:, 3) == 0;
d1 = data(l, 1 : 2)';
d2 = data(~l, 1 : 2)';
plot(d1(1, :), d1(2, :), 'r+');
plot(d2(1, :), d2(2, :), 'g+');
axis('equal');
 
p = 10; % kiek duomenų skirti apmokymui procentais
features = 2; % požymių kiekis
k = 3; % kiek kaimynų naudoti klasifikuojant
q = int16(size(data, 1) * p / 100); % kiek duomenų skirti apmokymui
training = data(1 : q, 1 : features); % apmokymo duomenys
plot(training(:, 1), training(:, 2), 'b.'); % pažymėti, kurie duomenys dalyvauja apmokyme
group = data(1 : q, features + 1); % apmokymo duomenų klasės
total = size(data, 1); % duomenų kiekis
correct1 = 0;
correct2 = 0;
for i = q + 1 : total % testuoti tik su tais duomenimis, kurie nenaudojami apmokymui
    g = data(i, features + 1);
    sample = data(i, 1 : features);
    c = knnclassify(sample, training, group, k);
    if c == g
        correct1 = correct1 + 1;
    else
        fprintf('Klaida: %d \n', i);
    end;
    c = myknnclassify(sample, training, group, k);
    if c == g
        correct2 = correct2 + 1;
    else
        fprintf('Klaida: %d \n', i);
    end;
    fprintf('Total iterations: %d, done: %d\n', total - q, i - q);
end;
fprintf('knnclassify rezultatas: %4.2f %%\n', 100 * correct1 / (total - q));
fprintf('myknnclassify rezultatas: %4.2f %%\n', 100 * correct2 / (total - q));

KNN funkcija:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

function [class] = myknnclassify(sample, training, group, k)
A = zeros(1, 2); % matrica artimiausiems kaimynams saugoti
n = size(training, 1); % iš viso kaimynų kiekis
for i = 1 : n % kiekvienam kaimynui atliekami skaičiavimai
    d = euclidean_distance(training(i, :), sample); % randamas Euklido atstumas
    s = size(A, 1);
    if s < k % jeigu artimiausių kaimynų atrinkta mažiau nei K, tai papildom
        A(s + 1, 1 : 2) = [group(i) d];
    elseif A(k, 2) > d % jeigu artimiausių kaimynų atrinkta K, tolimiausią keičiam artimesniu
        A(k, 1 : 2) = [group(i) d];
    end;
    A = sortrows(A, 2); % išrūšiuojam pagal atstumą
end;
m = max(group) + 1; % klasių kiekis
B = zeros(1, m);
for i = 1 : k % randam po kiek artimiausių kaimynų priklauso kiekvienai klasei
    B(A(i, 1) + 1) = B(A(i, 1) + 1) + 1;
end;
[maxC, maxI] = max(B); % klasės indeksas, kurioje daugiausiai artimiausių kaimynų
[minC, minI] = min(B); % klasės indeksas, kurioje mažiausiai artimiausių kaimynų
if maxC > minC % laimėtoja klasė, kuri turi daugiausia artimiausių kaimynų
    class = maxI - 1;
else % jeigu visos klasės turi vienodą kaimynų kiekį (gali būti ir kitokie kriterijai)
    class = -1; % atmetimo (reject) klasė
end;

Euklido atstumas:

1
2
3
4
5
6
7

function [d] = euclidean_distance(x, y)
n = size(x, 2);
d = 0;
for i = 1 : n
    d = d + (x(i) - y(i)) ^ 2;
end;
d = d ^ 0.5;

Abejų funkcijų rezultatai yra gana panašūs. Kadangi duomenų klasės nepersidengia ir yra gana toli viena nuo kitos, tai pakanka vos kelių procentų apmokymo duomenų, kad visus likusius duomenis abi funkcijos klasifikuotų šimto procentų tikslumu. Žemiau esančiame paveiksliuke matosi mėlynai pažymėti taškai, kurie buvo naudojami apmokymui ir yra atrinkti atsitiktine tvarka (vos 10 procentų nuo visų duomenų). Kiekvieną kartą užkrovus duomenis, parenkama vis kita 10 procentų duomenų imtis apmokymui, taip mes galime matyti kaip funkcijų tikslumas priklauso nuo skirtingų duomenų aibių, tačiau išlaikant norimą jos dydį:

myknnclassify funkcija yra dešimt kartų lėtesnė nei Matlab standartinė funkcija knnclassify ir kiek testavau, mažiau tiksli. Be to, aprašytoje gali būti (ir tikriausiai yra) klaidų, nes ji nėra pilnai ištestuota kaip Matlab funkcija. Todėl Matlab funkcija yra gerokai pranašesnė, kuri be to dar turi ir įvairių papildomų parametrų.

Tyrinėjimai artimiausių kaimynų metodas, euklido atstumas, k nearest neighbors, klasifikavimas, KNN

Apie pačią matematinę morfologiją plačiai nerašysiu, mes tik pabandysime pažiūrėti kelias sąvokas ir iš karto šį bei tą išbandyti.

Struktūrizavimo elementas (structuring elements) yra binarinis paveiksliukas, pavyzdžiui kokia nors geometrinė forma (stačiakampis, diskas, kaimynų aibė ir t.t.), dažniausiai gerokai mažesnis už paveiksliuką kuriam bus naudojamas.

Morfologiniai operatoriai:

• Išplėtimas (dilation) – padidina paveiksliuką;
• Erozija (erosion) – sumažina paveiksliuką;
• Uždarymas (closing) – uždaro vidines skyles ar išorines įlankas;
• Atidarymas (opening) – pašalina smulkias, atsikišusias į foną dalis.

Paimkim tokį pavyzdį:

Tai binarinis 10×10 dydžio paveiksliukas su figūra.

Paimkime paprastą struktūrizavimo elementą 3×3 dydžio (kvadratą): [1 1 1; 1 1 1; 1 1 1] ir panaudokime visus keturis operatorius (viršutinis kairysis – dilation, viršutinis dešinysis – closing, apatinis kairysis – erosion, apatinis kairysis – opening):

Kodas:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

tic;
clc;
close all;
file = 'D:\shape.bmp';
rgb = imread(file);
B = im2bw(rgb, 0.8);
 
se = strel('square', 3);
B1 = imdilate(B, se);
subplot(2,2,1), subimage(B1)
 
se = strel('square', 3);
B2 = imerode(B, se);
subplot(2,2,2), subimage(B2)
 
se = strel('square', 3);
B3 = imclose(B, se);
subplot(2,2,3), subimage(B3)
 
se = strel('square', 3);
B4 = imopen(B, se);
subplot(2,2,4), subimage(B4)
 
fprintf('Užtruko: %4.2f s\n', toc);

Naudojant šias operacijas atliekamos įvairiausios užduotys, pavyzdžiui galime išryškinti šių dvejų detalių defektus:

Jei labai gerai įsižiūrėti, tai abi detalės turi po defektą t.y. trūksta po vieną dantuką, tačiau tam kad pastebėti tokius defektus, reikia labai daug atidumo. Pabandysime naudojant morfologinius operatorius, šiuos defektus išryškinti. Kodas:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

tic;
clc;
close all;
file = 'D:\gears.bmp';
rgb = imread(file);
B = im2bw(rgb, 0.8);
imshow(B); 
 
% rasti skyles
file = 'D:\hole_ring.bmp';
rgb = imread(file);
hole_ring = im2bw(rgb, 0.8);
SE = strel(hole_ring);
B1 = imerode(B, SE); % erozija naudojant žiedo struktūros kaimynus
imshow(B1); 
 
% padidinti skyles
SE = strel('disk', 20, 0);
B2 = imdilate(B1, SE);
imshow(B2); 
 
% užkimšti skyles
B3 = bitor(B, B2);
imshow(B3); 
 
% tik dantukai
SE = strel('disk', 20, 0);
B7 = imopen(B3, SE);
B8 = B3 - B7;
imshow(B8);
 
% padidinti dantukai
SE = strel([1; 1]);
B9 = imerode(B8, SE);
SE = strel('disk', 4, 0);
B9 = imdilate(B9, SE);
imshow(B9);

Rezultatai:

Paskutiniame paveiksliuke defektai išryškinti akivaizdžiai. Jei žiūrėti paeiliui, tai pirmiausia erozijos pagalba aptinkamos skylės. Tam pritaikytas erozijos operatorius ir žiedo kaimynai, kurie pateikti atskiru paveiksliuku – hole_ring.bmp:

Gauti keturi maži objektai (antras paveiksliukas), kurie yra skylių centruose. Jie padidinami tiek, kad pilnai uždengtų pačias skyles. Naudojant OR operaciją su originaliu paveiksliuku ir tais gautais objektais, gaunamas trečias paveiksliukas be skylių. Ketvirtame atmesta vidinė detalės dalis, o penktame kiekvienas dantukas išryškinamas.

Tyrinėjimai atidarymas, closing, dilation, erosion, erozija, išplėtimas, matematiniai morfologiniai operatoriai, opening, uždarymas

Objektų žymėjimas (connected components labeling, objects labeling) yra atliekamas visus objektus vieną nuo kito išskiriant, pavyzdžiui nudažant skirtinga spalva arba priklijuojant etiketę (label). Išbandysim abu šiuos variantus.

Paimkim binarinį paveiksliuką:

Algoritmas (recursive labeling):

• skenuojamas visas paveiksliukas ieškant juodos spalvos ir jeigu tokia rasta – nuspalvosime išskirtine intensyvumo spalva (pilku atspalviu);
• sudaromas rasto taško kaimynų sąrašas ir pas kiekvieną jo kaimyną apsilankoma, nudažant juos visus vienodu intensyvumu. Lankyti kaimynus galima su rekursija, tuomet visi objekto taškai bus tikrai apeiti ir pažymėti;
• ieškomas naujas objektas.

Programos kodas žemiau:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

tic;
clc;
close all;
file = 'D:\objektai.bmp';
rgb = imread(file);
image(rgb);
global I;
I = rgb2gray(rgb);
label = 0;
labelDelta = 20;
[MaxRow, MaxCol] = size(I);
% set(0,'RecursionLimit',1500); didesniems objektams gali reikėti didinti
% rekursijų kiekio limitą
for L = 1 : MaxRow
    for P = 1 : MaxCol
        if I(L, P) == 0 % aptiktas objektas
            label = label + labelDelta; % kiekvienam objektui kitoks intensyvumas
            search(label, L, P); % nuspalvoti visus objekto pikselius vienodu intensyvumu
        end;
    end;
end;
imshow(I);
disp(toc);

Rekursija apeinanti visus kaimynus:

1
2
3
4
5
6
7
8
9
10
11
12

function [] = search(label, L, P)
global I;
I(L, P) = label; % keičiamas pikselio intensyvumas
[Nset] = neighbors(I, L, P); % randami 4 pikselio kaimynai
s = size(Nset);
n = s(3);
for i = 1 : n % kiekvienas kaimynas tikrinamas atskirai
    value = Nset(:, :, i);
    if value(3) == 0
        search(label, value(1), value(2));
    end;
end;

Galima rasti ir 8 ir 4 kaimynus. Šiuo atveju randami tik keturi:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

function [res] = neighbors(I, L, P)
[MaxRow, MaxCol] = size(I);
res = [];
cur = -1;
if L > 1 % viršutinis kaimynas
    cur = cur + 1;
    res(:, :, cur + 1) = int16([L - 1 P int16(I(L - 1, P))]);
end;
if L < MaxRow % apatinis kaimynas
    cur = cur + 1;
    res(:, :, cur + 1) = int16([L + 1 P int16(I(L + 1, P))]);
end;
if P > 1 % kairysis kaimynas
    cur = cur + 1;
    res(:, :, cur + 1) = int16([L P - 1 int16(I(L, P - 1))]);
end;
if P < MaxCol % dešinysis kaimynas
    cur = cur + 1;
    res(:, :, cur + 1) = int16([L P + 1 int16(I(L, P + 1))]);
end;

Apdorojus pirmąjį paveiksliuką:Algoritmas bėga per visus pikselius nuo viršaus apačion ir nuo kairės dešinėn, todėl pirmiausia aptinkamas trečias objektas ir po to visi kiti. Rezultate matosi kaip kiekvieno objekto intensyvumas skiriasi.

Pabandom lygiai tą patį padaryti su paveiksliuku, kurį naudojom skaičiuojant objektus, tik intensyvumo pokytį pakeičiam į 4, nes objektų gerokai daugiau ir reiktų, kad visi sutilptų į 255 reikšmes:

Labiausiai skiriasi vienas nuo kito objektai, kurie buvo pažymėti pirmiausiai ir paskiausiai. Tokį paveikslą gali būti sunku naudoti, nes nėra kaip pasakyti apie kurias ląsteles norima diskutuoti, juk nesakysi, kad „ląstelė, kurios intensyvumas 74 % …“, todėl geriausia, jei kiekvienam objektui suteiksim identifikatorius. Kodas žemiau (funkcijos nekeistos):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

tic;
clc;
close all;
file = 'D:\FROG BLOOD binary.bmp';
rgb = imread(file);
file2 = 'D:\FROG BLOOD 402 rgb.bmp';
rgb2 = imread(file2);
imshow(rgb2);
global I;
I = rgb2gray(rgb);
label = 0;
labelDelta = 4;
[MaxRow, MaxCol] = size(I);
% set(0,'RecursionLimit',1500); didesniems objektams gali reikėti didinti
% rekursijų kiekio limitą
for L = 1 : MaxRow
    for P = 1 : MaxCol
        if I(L, P) == 0 % aptiktas objektas
            label = label + labelDelta; % kiekvienam objektui kitoks intensyvumas
            search(label, L, P); % nuspalvoti visus objekto pikselius vienodu intensyvumu
            text('units','pixels','position',[P MaxRow-L],'fontsize',7,'string',num2str(label/labelDelta),'BackgroundColor',[0.9 0.9 0.9])
        end;
    end;
end;
% imshow(I);
disp(toc);

Rezultatas:

Visi 59 objektai pažymėti.

Šio algoritmo didžiausias trūkumas, kad netinka dideliems objektams, nes tada reikia didinti rekursijų kiekio limitą. O ir padidinus ne visada to pakanka, nes viršijus Matlab galimybes, jis lūžta. Tačiau šiai problemai spręsti yra ir daugiau algoritmų: “row-by-row labeling”, naudojant “union-find structure” arba “run-length encoding” ir kt.

Tyrinėjimai connected components labeling, objects labeling, objektų išskyrimas, objektų pažymėjimas, paveiksliukų apdorojimas