Kirjautuminen

Haku

Tehtävät

Keskustelu: Ohjelmointikysymykset: C++: Tehokas lokaalien maksimien ja minimien etsiminen

Sivun loppuun

msdos464 [22.03.2012 23:11:56]

#

Minulla on 4 "kerrosta" 1280 x 720 float arvoja, ja haluan löytää niiden lokaalit minimit ja maksimit mahdollisimman nopeasti Intel i5-2500 prosessorilla. Pisteen vertailu tehdään 8 naapuriin samassa tasossa ja 5 naapuriin "edellisessä" ja "seuraavassa" tasossa, eli yhteensä on enintään 8 + 2 * 5 = 18 vertailua. Tämä etsintä tehdään luonnollisesti kerroksille numerot 2 ja 3.

Oikeastaan lisäehtona että paikallisen minimin pitää olla vähintään tietyn muuttuvan toleranssin "eps" verran pienempi kuin sitä ympäröivien lukujen. Vastaavasti maksimin pitää olla "eps" verran isompi.

Data on float* taulukossa, jossa luku kohdassa (i, j) on löydyy paikasta data[j * width + i]. Nykyinen toteutus käy dataa läpi neljällä threadilla, joista kaksi käsittelee kerrosta 2 ja kaksi muuta käsittelee kerrosta 3. Ensimmäinen threadi käsittelee 2. kerroksen ylempää puoliskoa ja toinen threadi käsittelee sen alempaa puoliskoa. Iterointi tehdään riveittäin ja sarakkeittain järjestyksessä.

Kaikki algoritmit ovat O(n*m), mutta mietin että millainen access pattern olisi L1 ja L2 cachejen kannalta paras. Ehkä kannattaisi iteroida "sub-blockeja" kerrallaan? Pitää kokeilla huomenna käytännössä. Tietorakenteista on tuskin hyötyä, koska sen rakentamiseen kuluisi varmaan vähintään O(n*m*log(n*m)) aikaa.

Data on harmaasävykuvan laplacian-of-gaussian konvoluutiosta, eli paikallisia minimejä ja maksimeja on yleensä muutama sata.

Lisäys:

Tässä nykyinen koodi minimin tarkastamiseksi:

const float currentMin = current[j][i] + grayMargin;

if (
		currentMin < current[j][i-1]	&&
		currentMin < current[j][i+1]	&&
		currentMin < current[j-1][i]	&&
		currentMin < current[j+1][i]	&&
		currentMin < prev[j][i]			&&
		currentMin < next[j][i]			&&
		currentMin < current[j-1][i-1]	&&
		currentMin < current[j-1][i+1]	&&
		currentMin < current[j+1][i-1]	&&
		currentMin < current[j+1][i+1]	&&
		currentMin < prev[j+1][i]		&&
		currentMin < prev[j-1][i]		&&
		currentMin < prev[j][i+1]		&&
		currentMin < prev[j][i-1]		&&
		currentMin < next[j+1][i]		&&
		currentMin < next[j-1][i]		&&
		currentMin < next[j][i+1]		&&
		currentMin < next[j][i-1]) {
	// paikallinen minimi
}

Tuo vertailujärjestys tuntuu olevan tärkeä, jos esimerkiksi vertaillaan ensiksi "prev" ja "next" naapureihin, nopeus laskee muistaakseni 25%.

tkok [22.03.2012 23:55:13]

#

Suuri määrä rinnakkain suoritettavia float-operaatioita kannattaa suorittaa vektorirekistereillä vektorisuorittimessa. En ole varma löytyykö sieltä myös vertailuja suorittavia osia.

http://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html sanoo että:

Vectors are compared element-wise producing 0 when comparison is false and -1 (constant of the appropriate type where all bits are set) otherwise. Consider the following example.

     typedef int v4si __attribute__ ((vector_size (16)));

     v4si a = {1,2,3,4};
     v4si b = {3,2,1,4};
     v4si c;

     c = a >  b;     /* The result would be {0, 0,-1, 0}  */
     c = a == b;     /* The result would be {0,-1, 0,-1}  */

Eli ilmeisesti kyllä.
SSE:n käytöstä lisää: http://www.cs.uaf.edu/2009/fall/cs301/lecture/11_13_sse_intrinsics.html
Suorittimestasi löytyy AVX (Advanced Vector Extensions), jossa vektori rekisterit on kasvatettu 128bit->256bit joten sinne pitäisi mahtua 8 floattia(32bit) rinnan kerralla, http://en.wikipedia.org/wiki/Advanced_Vector_Extensions
EDIT: Lisäksi tämä kannattaa katsoa: http://www.agner.org/optimize/optimizing_cpp.pdf Luku 12 käsittelee vektoreita ja muutenkin hyvä paperi optimoinnista.

msdos464 kirjoitti:

Tuo vertailujärjestys tuntuu olevan tärkeä, jos esimerkiksi vertaillaan ensiksi "prev" ja "next" naapureihin, nopeus laskee muistaakseni 25%.

Mitä kääntäjää käytät ja millä optimointiasetuksilla? Tuo vaikuttaisi siltä, että kääntäjäsi ei tunnista tuossa mahdollisuutta käyttää vektorointia ja jättää sen optimoimatta. Toinen mahdollisuus on myös, että ehtolauseesi osoittautuu sattumalta nopeammin epätodeksi eri järjestyksellä.

msdos464 [23.03.2012 00:19:46]

#

Käyttis on 64 bittinen Ubuntu 10.04, g++ -O3 asetuksilla, en muista ulkoa g++ versiota mutta se mahtaa olla tuorein. Tiesin SSE:n olemassaolosta mutta en jotenkin tajunnut että sillä voisi olla mahdollista tehdä myös vertailuoperattori.

Luulen nopeuseron syynä olevan se, että current[j][ i-1 ] ja current[j][ i+1 ] ovat todennäköisemmin cacheissa kuin nuo prev ja next rivien arvot.

Täytyypä lueskella noista lisää. Tuo optimizing_cpp.pdf vaikuttaa myös hyvältä.

Tässä esimerkki dataa jos jotakuta kiinnostaa, 1280 x (4 * 720) kuva tyypillisestä datasta. Ylin on alkuperäinen kuva, ja sitä seuraavat eriasteisia LoG suodatuksia.

Metabolix [23.03.2012 00:56:52]

#

Olisikohan kuitenkin myös algoritmissa optimoinnin varaa? Yksinkertaisena ideana esitän, että koska yhdessä neljän ruudun neliössä voi olla vain yksi paikallinen minimi, voisit vähentää silmukan iteraatiot neljännekseen ja ennen vertailua valita näistä neljästä tutkittavaksi vain parhaan. Ainakin omassa testissäni, jossa täytin taulukon satunnaisdatalla, tämä ratkaisu nopeutti ohjelmaa 25 %.

Vastaavasti jos suurin osa taulukosta koostuu pitkistä monotonisista osuuksista, voisi auttaa tiivis silmukka, joka rullaisi kyseisen monotonisen osuuden loppuun, siis while (taulu(i+1) < taulu(i)) ++i. Satunnaisdatallani tämä tietenkin vain hidasti.

Onko tuossa mallidatassasi siis tarkoitus vertailla kuvia juuri tuossa järjestyksessä? Selittäisitkö vielä, mitä tällä saavutetaan? :)

Taas kerran clang++ tuottaa paremman tuloksen kuin g++.

msdos464 [23.03.2012 01:05:46]

#

Kyseessä on Scale-invariant Feature Transform", jolla löydetään kuvista stabiileja "piirrepisteitä".

Tällä on lukuisia sovelluskohteita kuten objektien tunnistus, mutta itse käytän sitä videodataan. Tuolla voi toteuttaa vaikka digitaalisen kuvanvakauksen, mutta itse aion paikantaa kameran paikkaa maailmassa. Tämä tulee olemaan diplomityöni.

Tässä on 60 mb video vanhahkosta versiosta, joka laski "optista vuota" 8 fps. Nykyinen yltää jo 20 fps nopeuteen paremmalla tarkkuudella.

Täällä on esimerkki tunnettujen objektien etsimisestä valokuvista.

Opas sanoo "Don't change && to & unless you expect the && expression to
generate many branch mispredictions.", hmmm tuossahan noita on iso läjä, ja branch factor on melkoinen. Voisi koittaa muuttaa osan noista operaatioista bitti-andeiksi, ja ryhmitellä siten että cachesta olisi eniten hyötyä.

Ooh muistin että jokainen noista current[j][i] on funktiokutsu tähän:

const float *operator[](const int rowIndx) const {
	return _destination + rowIndx * _w;
}

Kannattaisi varmaan tallentaa nuo rivien alkukohdat talteen, koska en yhtään tiedä että miten kääntäjä on optimoinut tuon :D

FooBat [23.03.2012 05:08:02]

#

Yritä pitää mahdollisimman paljon luettua dataa rekistereissä sen sijaan, että jokaisessa vertailussa luet datat uudestaan muistista. Jos vertailet silmukassa vuoron perään vieräkkäisiä pisteitä, tulet lukeneeksi samat datat moneen kertaan.

Sijoita luetut arvot lokaaleihin muuttujiin (yleensä pysyvät rekistereissä) ja kierrätä niitä sopivasti siten, että vain uudet reunapikselit pitää lukea muistista.

msdos464 [23.03.2012 09:50:55]

#

Hmm outoa, jos tallennan noiden rivien alkupointterit (yhteensä 3 kpl) "const float*" muuttujiin, koodi hidastuu 25%. Ehkä ne sitten vievät tilaa rekistereistä.

Äh en tiedä johtuiko se siitä... Nyt kun kikkailin kääntäjän parametreilla ja palasin alkuperäiseen koodiin niin sekin näyttää hidastuneen.

Muok: Ei kun se nopeutui 25%... luin nuo tuolokset ristiin. -.-

Epsilon [23.03.2012 21:03:58]

#

Kokeilin nopeasti SSE-käskykannan kanssa määrittää minimiä neljälle vierekkäiselle pikselille kerrallaan. Valitettavasti suorittimestani ei löydy AVX-käskykantaa, joten en voinut testata sillä. (Pitäisi kai päivittää vähän uudempaan rautaan..) Kokeilusoftan käänsin MS:n Windows SDK 7:n mukana tulevalla 64-bittisellä compilerilla, joten en tiedä, millä kaikilla kääntäjillä nuo intrinsicit toimivat sellaisenaan.

Kokeilussa alla oleva koodinpätkä tuotti yhdessä threadissa ajettuna noin 2.8x speedupin alkuperäisessä viestissä mainittuun "if(monta comparea and:lla yhdistettynä)" -versioon verrattuna. Tuskin tuo kokeilemani versio on mitenkään optimaalinen, mutta varmaan tämän rohkaisemana kannattaa ainakin kokeilla niitä AVX-käskyjä ja katsoa millaista nopeusparannusta niiden kanssa on saatavissa aikaan.

(Toivottavasti en typottanut mitään kovin pahasti. Kokeilutuloksen jatkokäyttö kuitenkin joka tapauksessa omalla vastuulla.. =P)

Edit: Huom! Noissa loopeissa ei ole mukana sitä grayMargin -thresholdin lisäämistä, mutta se ei nyt kai ollut oleellinen pointti kokeilun kannalta.

const int WIDTH  = 1280;
const int HEIGHT = 720;

inline int c2i(int i, int j) { return j*WIDTH + i; }

void FindLocalMinima_sse(const float * prev, float * current, float * next, int * result)
{
    const int H = HEIGHT - 1;
    const int W = WIDTH - 1;
    for(int j = 1; j < H; ++j)
    {
        int i = 1;
        for(; i < W-3; i+=4)
        {
             __m128 xmm_current_min = _mm_loadu_ps(current + c2i(i, j));

             __m128 xmm_result  =                        _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i-1, j  )));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i+1, j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i  , j-1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i  , j+1))));

             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i-1, j-1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i+1, j-1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i-1, j+1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(current + c2i(i+1, j+1))));

             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(next + c2i(i  , j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(next + c2i(i-1, j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(next + c2i(i+1, j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(next + c2i(i  , j-1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(next + c2i(i  , j+1))));

             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(prev + c2i(i  , j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(prev + c2i(i-1, j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(prev + c2i(i+1, j  ))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(prev + c2i(i  , j-1))));
             xmm_result         = _mm_and_ps(xmm_result, _mm_cmplt_ps(xmm_current_min, _mm_loadu_ps(prev + c2i(i  , j+1))));

             _mm_storeu_ps((float *) result + c2i(i, j), xmm_result);
        }
        for(; i < W; ++i)
        {
            const float currentMin = current[c2i(i, j)];

            if (currentMin < current[c2i(i-1, j  )]	&&
                currentMin < current[c2i(i+1, j  )]	&&
                currentMin < current[c2i(i  , j-1)]	&&
                currentMin < current[c2i(i  , j+1)]	&&
		currentMin < prev[c2i(i, j)]    &&
                currentMin < next[c2i(i, j)]	&&
                currentMin < current[c2i(i-1, j-1)]	&&
                currentMin < current[c2i(i+1, j-1)]	&&
                currentMin < current[c2i(i-1, j+1)]	&&
                currentMin < current[c2i(i+1, j+1)]	&&
                currentMin < prev[c2i(i  , j-1)]	&&
                currentMin < prev[c2i(i  , j+1)]	&&
                currentMin < prev[c2i(i-1, j  )]	&&
                currentMin < prev[c2i(i+1, j  )]	&&
                currentMin < next[c2i(i  , j-1)]	&&
                currentMin < next[c2i(i  , j+1)]	&&
                currentMin < next[c2i(i-1, j  )]	&&
                currentMin < next[c2i(i+1, j  )] )
            {
	            result[c2i(i, j)] = -1;
            }
            else
            {
                result[c2i(i, j)] = 0;
            }
        }

    }
}

tepokas [23.03.2012 22:25:09]

#

msdos464 kirjoitti:

hmmm tuossahan noita on iso läjä, ja branch factor on melkoinen

kokeile muotoa

t[0]=current[][]
t[1]=...
t[n]=next[][]

j=0;
for(ii=0 to n) j=j+(0<t[ii]-currentMin);

if(j/(0<1)==(n+1)) -> minimi

ja voit saada poistettua branchejä

Muok. BugiFix 0.1
Mod. korjasi kooditagit!

User137 [23.03.2012 22:32:39]

#

Sitten voi olla myös mahdollista käyttää pienempää hakuresoluutiota, esim skippaa joka toisen pikselin. Tai jos löytää minimiarvon niin tutkii sitten vierestä skipattukin pikseli, josko siinä olisi parempi arvo.

Epsilon [24.03.2012 10:04:46]

#

Tarkentava kysymys: Alkuperäisessä postissa puhut vain neljän floatin kerroksesta, mutta onko oikeasti tilanne se, että aiot tehdä analyysia oikeastaan videostreamille? Eli siis frameja olisi tarjolla oikeastaan jatkuva virta eikä toisistaan irrallisia (prev, current, next) -kolmikoita? Videostreamilla tarkoitan tässä siis sitä, että jokaisen prosessointikierroksen alussa asetetaan prev=current, current=next ja next=(uusi data), eli käytännössä aina mukaan tulee yksi uusi frame ja aiempia käytetään sellaisenaan. (Tuosta linkkaamastasi videosta tuli mieleen, että saattaisit tehdä jotain tämän tyylistä..)

Videostreamia ajatellen tulee mieleen, että voisi yrittää cachettaa jokaisen framen jokaiselle pikselille preprosessointivaiheessa (vaikka silloin kun filtteröit sitä kuvaa) naapuruston {(-1, 0), (1, 0), (0, -1), (0, 1)} minimiarvon. Tällöin varsinaista minimiä etsittäessä saisit korvattua 12 vertailua kolmella vertailulla, eli alkuperäinen iffisi muuttuisi suunnilleen alla olevaan muotoon:

if( //Vertailu valmiiksi laskettujen naapurustojen minimeihin
    currentMin < prev_ngbr_min[c2i(i, j)] &&
    currentMin < current_ngbr_min[c2i(i, j)] &&
    currentMin < next_ngbr_min[c2i(i, j)] &&

    //prev ja next -naapurustojen keskipisteet
    currentMin < prev[c2i(i, j)]  &&
    currentMin < next[c2i(i, j)]  &&

    //Diagonaalin suunnassa olevat "current"-naapuripikselit,
    //jotka eivät ole mukana minimicachetuksessa
    currentMin < current[c2i(i-1, j-1)] &&
    currentMin < current[c2i(i+1, j-1)] &&
    currentMin < current[c2i(i-1, j+1)] &&
    currentMin < current[c2i(i+1, j+1)]
    )

Ja tuo iffi tietysti rinnakkaistuu vastaavasti SSE:llä (tai AVX:llä) kuin aiempikin. Kokeilin nopeasti aiemmin esittämääni SSE-looppiin tuon tyylistä viritystä ja käytännössä vertailuloopista tuli noin hieman vajaa 2x nopeampi yllä mainituilla virityksillä. Tämä ei tietysti yllätä, koska kaikkia operaatioita, muistihaut mukaan lukien tarvitaan vain noin puolet suhteessa alkuperäiseen. Tuon toteuttamalla liikuttaisiin SSE:llä siis jossain 5.5-6x nopeudessa esittämääsi alkuperäiseen if(...)-versioon nähden. Tässä nopeutuksessa ei tosin ole huomioitu sitä, että preprosessoinnissa joudutaan jokaiselle framelle tekemään vähän raskaampi operaatio, kun etsitään nuo naapuruston minimit.

Tietysti tuota minimien esimääritystä voi jatkaa vielä siten, että tallettaa jokaiselle pikselille esimerkiksi naapurustojen {(-1, 0), (1, 0)} miniarvon vielä yhteen taulukkoon. Tuo minimi tulee efektiivisesti laskettua jo silloin, kun lasketaan noita koko naapuruston minimejä (tämä on vain osa siitä), joten varsinaista lisävaivaa laskennasta ei tule. Tällöin tätä naapurustotaulukkoa käyttämällä voi eliminoida vielä kaksi noista iffin sisällä olevista compareista, eli lopputulosiffi olisi jotain tyyliin:

if( //Vertailu valmiiksi laskettujen naapurustojen minimeihin
    currentMin < prev_ngbr_min[c2i(i, j)] &&
    currentMin < current_ngbr_min[c2i(i, j)] &&
    currentMin < next_ngbr_min[c2i(i, j)] &&

    //prev ja next -naapurustojen keskipisteet
    currentMin < prev[c2i(i, j)]  &&
    currentMin < next[c2i(i, j)]  &&

    //Diagonaalin suunnassa olevat "current"-naapuripikselit,
    //jotka eivät ole mukana minimicachetuksessa
    currentMin < current_ngbr2_min[c2i(i, j-1)] &&
    currentMin < current_ngbr2_min[c2i(i, j-2)]
    )

Esimerkiksi SSE:n kannalta tuo iffi näyttäisi jo sillä tavalla kivalta, että preprosessoinnissa voisi järjestää kaikkien taulukoiden aligmentin siten, että nuo rivin i:nnet alkiot (i:n incrementti aina 4) osuisivat aina 16-byte aligned kohtaan, jolloin SSE:llä ei tarvitsisi ladata unaligned muistipaikasta ensin tavaraa rekisteriin vaan SSE-operaatiot (esim. alla olevan esimerkin min-operaatiot) pystyisivät lukemaan tavaran suoraan 16-byte aligned muistipaikasta. Tuosta aligment-asian huomioinnista saisi varmaan vielä jonkunmoista boostia tuohon varsinaiseen innerlooppiin.

Ja vielä kevyt "kasvojen kohotus" tuohon aiemmin lähettämääni SSE-viritykseen. Käytin siinä turhaan and:ia ja comparea joka välissä. Vähemmillä käskyillä pääsee, kun käyttää vain min-operaatiota toistuvasti ja tekee vain yhden comparen lopussa. Tuolla virityksellä tuli vajaa 10% lisää nopeutta. Parannettu versio alla (vain sse-osuus loopin sisältä, parannettu vain tuon compare aspektin suhteen, ei tehty em. preprosessoinnin hyödyntämistä):

...
__m128 xmm_current_min = _mm_loadu_ps(current + c2i(i, j));
__m128 xmm_ngbr_min  =                           _mm_loadu_ps(current + c2i(i-1, j  ));
__m128 xmm_ngbr_min2 =                           _mm_loadu_ps(current + c2i(i+1, j  ));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,  _mm_loadu_ps(current + c2i(i  , j-1)));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2, _mm_loadu_ps(current + c2i(i-1, j-1)));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,  _mm_loadu_ps(current + c2i(i+1, j-1)));

xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(current + c2i(i  , j+1)));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(current + c2i(i-1, j+1)));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(current + c2i(i+1, j+1)));

xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(next + c2i(i  , j  )));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(next + c2i(i-1, j  )));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(next + c2i(i+1, j  )));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(next + c2i(i  , j-1)));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(next + c2i(i  , j+1)));

xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(prev + c2i(i  , j  )));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(prev + c2i(i-1, j  )));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(prev + c2i(i+1, j  )));
xmm_ngbr_min         = _mm_min_ps(xmm_ngbr_min,   _mm_loadu_ps(prev + c2i(i  , j-1)));
xmm_ngbr_min2        = _mm_min_ps(xmm_ngbr_min2,  _mm_loadu_ps(prev + c2i(i  , j+1)));

_mm_storeu_ps((float *) result + c2i(i, j), _mm_cmplt_ps(xmm_current_min, _mm_min_ps(xmm_ngbr_min, xmm_ngbr_min2)));
...

Edit: Tarkentava kysymys taisi olla vähän hölmö, kun luin alkuperäisen viestiketjun tarkemmin ja löysin jopa ihan kuvan, jossa nuo mainitut neljä kerrosta oli näytetty. Jätän kuitenkin tämän pohdinnan tähän, jos vaikka siitä tulisi mieleen joitain käyttökelpoisia ajatuksia..

FooBat [24.03.2012 14:22:26]

#

Alkuperäiseen cache kysymykseen vastaus on, että kannattaa tehdä kaikki lokaalit vertailut peräkkäin. Olettaisin, että seuraava on aika hyvä. Seuraavan ja edellisen kuvan indeksöintiä kannattaa välttää mahdollisimman paljon, koska todennäköisesti niiden data ei ole cachessa. Cache miss on paljon suurempi penalty verrattuna branch missiin, joten ei kannata edes yrittää poistaa kaikkia branchejä tuosta.

float prevIten, currItem, nextItem;
prevItem = current[j][0];
currItem = current[j][0];


for i:

nextItem = current[j][i+1];
const float currentMin = currItem + grayMargin;

if (
		currentMin < prevItem	&&
		currentMin < nextItem	&&
		currentMin < current[j-1][i-1]	&&
		currentMin < current[j-1][i]	&&
		currentMin < current[j-1][i+1]	&&

		currentMin < current[j+1][i-1]	&&
		currentMin < current[j+1][i]	&&
		currentMin < current[j+1][i+1]	&&

		currentMin < prev[j][i-1]		&&
		currentMin < prev[j][i]			&&
		currentMin < prev[j][i+1]		&&

		currentMin < next[j][i-1]       &&
		currentMin < next[j][i]			&&
		currentMin < next[j][i+1]		&&



		currentMin < prev[j-1][i]		&&
		currentMin < prev[j+1][i]		&&
		currentMin < next[j-1][i]		&&
		currentMin < next[j+1][i]
		) {
	// paikallinen minimi
}
prevItem = currItem;
currItem = nextItem;

msdos464 [24.03.2012 15:07:46]

#

Äh en ehdi lukea viestejä kun pitää lähteä liikkeelle, mutta vastaan pikaisesti tähän:

Epsilon kirjoitti:

Alkuperäisessä postissa puhut vain neljän floatin kerroksesta, mutta onko oikeasti tilanne se, että aiot tehdä analyysia oikeastaan videostreamille?

Joo analyysi tehdään videostreamille, mutta nuo minimit ja maksimit etsitään ainoastaan yhdelle framelle kerrallaan. Myöhemmässä vaiheessa peräkkäisistä frameista etsitään toisiaan vastaavat piirrepisteet, jonka perusteella tuo "optinen vuo" voidaan laskea.

Vilkaise vielä tätä kuvaa: img_and_levels.png

Tuossa on ylhäällä alkuperäinen frame, ja sitten kolme eriasteista sumennusta. Tässä tapauksessa 3. kuvan (eli 2. sumennetun kuvan) pikseleitä vertailtaisiin ylimpään ja alimpaan sumennettuun kuvaan. Varsinaisessa sovelluksessa yhdestä kuvasta tehdään 4 eri sumennusversiota, koska käytössä on quadcore prosessori.

Luen muut viestit ajatuksella läpi myöhemmin. Kiitos ehdotuksista :)

Edit: Aikataulu muuttui vähän, niin ehdinkin lukea viestit.

Tässä on kuva jossa on alkuperäinen frame, ja siitä lasketut 4 "sumennettua" kuvaa: img_and_levels_2.png

Valitettavasti png on vain 8-bit harmaasävy, eli tuo ei ihan vastaa alkuperäistä dataa. Tässä on txt tiedosto jossa luvut ovat 8 desimaalin tarkkuudella: img_and_levels_2.zip

Tässä sovelluksessa on ikävää että data muuttuu noin 30 kertaa sekunnissa, eli mitään hienoja tietorakenteita ei monestikkaan ehdi rakentaa ja hyödyntää, verrattuna naiivimpaan lähestymistapaan.

Tuossa vaiheessa näytönohjain laskee CUDAlla edellisten kuvien piirteiden matchausta, joten sitä ei voida käyttää tässä vaiheessa algoritmia. Toisaalta näytönohjaimelle pitäisi kopioida melko suuri määrä dataa, jos nuo sumennukset lasketaan kuitenkin prossulla.

FooBat [24.03.2012 17:03:06]

#

Jos lasket nuo sumennukset itse aikaisemmassa vaiheessa, kannatta jo tähän vaiheeseen yhdistää esiprosessointia tätä toista vaihetta varten. Voit esim. jo esimmäisessä vaiheessa löytää yksittäisen kuvan lokaalit miniehdokkaat (minimi 8 ympäröivään nähden) ja tallentaa ne erilliseen listaan. Samoin jos et oikeasti tarvitse sumennuskuvia mihinkään muuhun kuin tähän vaiheeseen, voit tallentaa jokaiseen pikseliin viiden pikselin minimin, etkä oikeaa arvoa (saattaa vaati pientä välipuskuria, josta lasketaan minimit). Tällöin riitäisi käydä läpi vain listaan tallennettujen mahdolisten minimien kohdalta yksi tarkistus aikaisempaan ja seuraavaan tasoon.

Etuna tässä menettelyssä olisi, että pienellä lisätyöllä ensimmäisessä vaihessa voitaisiin pienentää toisen vaiheen muistinkäyttö ja cache missit minimiin. Ensimmäisessä sumennus vaiheessa joudut kuitenkin käsittelemään kaikkia kuvan pisteitä niin samalla kertaa ja samalla muisti accessilla voi tehdä myös hiukan lisätyötä.

Kannattaa myös miettiä tarvitaan dataa käsitellä 32 bittisinä floatteina. 16 intit puduttaisivat muistijäljen puoleen ja 8 bittiset neljäsosaan. Samalla cachen toimivuus paranisi suunnilleen samassa suhteessa.

msdos464 [24.03.2012 22:29:48]

#

Nykyinen koodi laskee tuon Laplacian of Gaussian sumennuksen (eli konvoluution) Fourier-tasossa, käyttäen FFTW kirjaston käänteis-FFT muunnosta. En siis käy pikseleitä yksi kerrallaan läpi tuota sumennusta laskettaessa.

Monet toteutukset käyttävät sellaista lähestymistapaa että laskevat ensiksi gaussisen sumennuksen (joka voidaan laskea rivi ja sarake kerrallaan) ja vähentävät kaksi tällaista sumennuksen tasoa toisistaan. Tämä on nimeltään Difference of Gaussian. Ehkä tuossa tapauksessa voitaisiin tehdä tällaista esiprosessointia tuota minimien ja maksimien etsintää varten.

Epsilon [24.03.2012 23:06:56]

#

Tein vielä vähän virittelyjä tuohon SSE-viritykseeni. Alla oleva koodi toimii samoin kuten edellinen versio, mutta se laskee kahta allekaista riviä kerrallaan. Koska allekaisilla riveillä on noissa naapurustoissa noin 1/3-osa yhteisiä minimiarvoa rajaavia pikseleitä, niin nämä tarvitsee hakea vain kerran.

void FindLocalMinima_sse2(const float * prev, float * current, float * next, int * result)
{
    const int H = HEIGHT - 1;
    const int W = WIDTH - 1;
    for(int j = 1; j < H-1; j+=2)
    {
        int i = 1;
        for(; i < W-3; i+=4)
        {
             //Luetaan current -kuvasta alue, jossa kaksi 3x3-naapurustoa leikkaavat
             __m128 xmm_common_min  =                            _mm_loadu_ps(current + c2i(i-1, j  ));
             __m128 xmm_row0_val    =                            _mm_loadu_ps(current + c2i(i  , j  ));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(current + c2i(i+1, j  )));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(current + c2i(i-1, j+1)));
             __m128 xmm_row1_val    =                            _mm_loadu_ps(current + c2i(i  , j+1));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(current + c2i(i+1, j+1)));

             //Luetaan current -kuvasta alue, joka kuuluu vain row0:n minimiin
             __m128 xmm_row0_min    = _mm_min_ps(xmm_row1_val,   _mm_loadu_ps(current + c2i(i-1, j-1)));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(current + c2i(i  , j-1)));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(current + c2i(i+1, j-1)));

             //Luetaan current -kuvasta alue, joka kuuluu vain row1:n minimiin
             __m128 xmm_row1_min    = _mm_min_ps(xmm_row0_val,   _mm_loadu_ps(current + c2i(i-1, j+2)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(current + c2i(i  , j+2)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(current + c2i(i+1, j+2)));

             //Luetaan prev-kuva (ristien keskimmäiset pikselit kuuluvat molempien rivien minimimaskiin)
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(prev + c2i(i  , j-1)));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(prev + c2i(i-1, j  )));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(prev + c2i(i  , j  )));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(prev + c2i(i+1, j  )));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(prev + c2i(i-1, j+1)));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(prev + c2i(i  , j+1)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(prev + c2i(i+1, j+1)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(prev + c2i(i  , j+2)));

             //Luetaan next-kuva (ristien keskimmäiset pikselit kuuluvat molempien rivien minimimaskiin)
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(next + c2i(i  , j-1)));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(next + c2i(i-1, j  )));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(next + c2i(i  , j  )));
             xmm_row0_min           = _mm_min_ps(xmm_row0_min,   _mm_loadu_ps(next + c2i(i+1, j  )));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(next + c2i(i-1, j+1)));
             xmm_common_min         = _mm_min_ps(xmm_common_min, _mm_loadu_ps(next + c2i(i  , j+1)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(next + c2i(i+1, j+1)));
             xmm_row1_min           = _mm_min_ps(xmm_row1_min,   _mm_loadu_ps(next + c2i(i  , j+2)));

             //Yhdistetään yhteiset minimit rivien erillisiin minimeihin
             xmm_row0_min           = _mm_min_ps(xmm_row0_min, xmm_common_min);
             xmm_row1_min           = _mm_min_ps(xmm_row1_min, xmm_common_min);

             //TODO: Lisää tässä xmm_row0_val:iin ja xmm_row1_val:iin se grayMargin threshold

             //Talleta tulokset
             _mm_storeu_ps((float *) result + c2i(i, j  ), _mm_cmplt_ps(xmm_row0_val, xmm_row0_min));
             _mm_storeu_ps((float *) result + c2i(i, j+1), _mm_cmplt_ps(xmm_row1_val, xmm_row1_min));

        }
        for(; i < W; ++i)
        {
            for(int j_inc = 0; j_inc < 2; ++j_inc)
            {
                const float currentMin = current[c2i(i, j+j_inc)];

                if (currentMin < current[c2i(i-1, j  +j_inc)]	&
                    currentMin < current[c2i(i+1, j  +j_inc)]	&
                    currentMin < current[c2i(i  , j-1+j_inc)]	&
                    currentMin < current[c2i(i  , j+1+j_inc)]	&

                    currentMin < prev[c2i(i, j+j_inc)]    &
                    currentMin < next[c2i(i, j+j_inc)]	  &

                    currentMin < current[c2i(i-1, j-1+j_inc)]	&
                    currentMin < current[c2i(i+1, j-1+j_inc)]	&
                    currentMin < current[c2i(i-1, j+1+j_inc)]	&
                    currentMin < current[c2i(i+1, j+1+j_inc)]	&

                    currentMin < prev[c2i(i  , j-1+j_inc)]	&
                    currentMin < prev[c2i(i  , j+1+j_inc)]	&
                    currentMin < prev[c2i(i-1, j  +j_inc)]	&
                    currentMin < prev[c2i(i+1, j  +j_inc)]	&

                    currentMin < next[c2i(i  , j-1+j_inc)]	&
                    currentMin < next[c2i(i  , j+1+j_inc)]	&
                    currentMin < next[c2i(i-1, j  +j_inc)]	&
                    currentMin < next[c2i(i+1, j  +j_inc)] )
                {
                    result[c2i(i, j+j_inc)] = -1;
                }
                else
                {
                    result[c2i(i, j+j_inc)] = 0;
                }
            }
        }

    }
}

Tuo sse-optimoitu koodinpätkä etsii koneellani (Core i7 920) 1280x720-kokoisille kuville minimit 336 kertaa sekunnissa, kun laskentaa suoritetaan yhdessä threadissa. Eli jos yhtä framea kohden noita samoja laskelmia pitäisi tehdä 2 kpl, niin silti ehdittäisiin käsitellä noin 170 framea sekunnissa. Eli noin 30 fps vauhdissa tuo minimien etsintä veisi karkeasti ottaen vähän yli kuudesosan framen ajasta.

msdos464 kirjoitti:

Nykyinen koodi laskee tuon Laplacian of Gaussian sumennuksen (eli konvoluution) Fourier-tasossa, käyttäen FFTW kirjaston käänteis-FFT muunnosta. En siis käy pikseleitä yksi kerrallaan läpi tuota sumennusta laskettaessa.

Monet toteutukset käyttävät sellaista lähestymistapaa että laskevat ensiksi gaussisen sumennuksen (joka voidaan laskea rivi ja sarake kerrallaan) ja vähentävät kaksi tällaista sumennuksen tasoa toisistaan. Tämä on nimeltään Difference of Gaussian. Ehkä tuossa tapauksessa voitaisiin tehdä tällaista esiprosessointia tuota minimien ja maksimien etsintää varten.

Veikkaisin, että tuo gaussian blurrin separoituvuutta kannattaisi tässä tapauksessa koittaa hyödyntää, jos kerran nopeus on kriittistä. FFT on käsittääkseni siinä suhteessa ikävä, että se generoi paljon randomia muistiaccessia kun taas separoituvilla filttereillä saa aikaan mukavaa lineaarista muistiaccessia. Lisäksi FFT:ssä pitää laskea kompleksiluvuilla, joten muistissa liikkuu äkkiä tuplamäärä tavaraa ja yksittäisen alkion laskutoimitukset ovat raskaampia.

Tietysti vähän filtterin koosta riippuu, missä vaiheessa "FFT+alkioittainen kertolasku" alkaa voittaa "separoitu + 1D-konvoluutio" toteutuksen, mutta ainakin noiden lähettämiesi kuvien perusteella näyttää siltä, että et käytä kauhean isolla radiuksella varustettuja filttereitä. Näin ollen "separoitu+konvoluutio"-implementaatio saattaisi käyttötilanteessasi olla FFT-versiota nopeampi. (En siis kokeillut asiaa vaan tämä on ainoastaan "mahanpohjafiilis" joka perustuu muutamiin aiempiin juttuihin, joissa olen tehnyt vähän vastaavan tyyppisiä filtteröintejä..)

msdos464 [24.03.2012 23:31:40]

#

Alunperin käytin FFT:tä koska aioin käyttää myös isompia sumennuskernelien kokoja. Nykyinen toteutus ei enää vaadi kauhean isoja kerneleitä, koska suuremmassa mittakaavassa ne voidaan korvata kuvan resoluutiota pienentämällä (niinkutsuttu kuvapyramidi).

Nykyisellä versiolla kuluu ajasta 25% eteenpäin muunnokseen (1 kpl), 50% ajasta neljään taaksepäinmuutokseen (neljällä säikeellä) ja 25% ajasta minimien ja maksimien etsimiseen. Kokonaisuudessaan framerate on noin 42.

Mikäli käyttäisin gaussiania (DoG) tuon LoGin sijaan, pitää huomata että tarvitsen yhden gaussisen sumennustason enemmän, jotta saisin tuloksena neljä difference of gaussian tasoa.

Täytynee tutustua noihin SSE ja AVX komentoihin.

Jaska [25.03.2012 11:12:10]

#

Nykyään FFT ei ole enää nopein tapa laskea muunnos. En osaa sanoa, paljonko paperissa http://arxiv.org/abs/1201.2501v1 esitetty algoritmi nopeuttaisi laskentaa.

msdos464 [25.03.2012 23:24:40]

#

Hmm, koska lasken tuossa kaistanpäästösuodatettuja signaaleja, käänteisfouriermuunnos voisi olla mahdollista toteuttaa nopeammassa ajassa kuin tuo FFTW toteutus. Täytyy katsoa että kuinka suuri osa noista komponenteista eroaa oleellisesti nollasta, että kannattaisiko tuota yrittää toteuttaa itse tehokkaammin.

Onko tuota uutta "Discrete Sparse Fourier Transform"ia saatavilla jo C kirjastona? Itse en onnistunut löytämään.

Jaska [26.03.2012 00:12:50]

#

msdos464 kirjoitti:

Onko tuota uutta "Discrete Sparse Fourier Transform"ia saatavilla jo C kirjastona? Itse en onnistunut löytämään.

Valitettavasti en tiedä. Luin algoritmista vain New Scientistista.


Sivun alkuun

Vastaus

Aihe on jo aika vanha, joten et voi enää vastata siihen.

Tietoa sivustosta