[C++] Vergelijken van afbeeldingen

[maleadt]

Ik ben momenteel een programma aan het schrijven waarvoor ik 2 afbeeldingen moet vergelijken, en een getal teruggeven die de orde van gelijkenis uitdrukt (hoger = sterker gelijkend). Onderdeel van de fitness() functie van een GA

Anyway, de afbeeldingen zijn Cairo surfaces, waarvan ik momenteel de raw data (unsigned char*) en die pixel voor pixel vergelijk. Elke pixel is 4 bytes in de reeks karakters (BGRA), waarvan ik voor elk kleur het verschil tussen beide kwadrateer, en per pixel die onder een wortel zet en bij een variabele optel (difference += sqrt( (red2-red1)² + (blue2-blue1)² + (green2-green1)²)). De uiteindelijke "resemblance" factor is dan 1 op die "difference" factor.
Ik ben hier op gekomen naar analogie van de gewone afstand tussen twee punten, maar ik weet niet in welke mate dit de beste methode is.

Nu heb ik al wat rondgezocht op internet, maar ik heb niet echt een algoritme gevonden dat echt gericht is op snelheid (of het moet een gewone memcpy geweest zijn, die te eenvoudig is in dit geval). Aangezien deze routine heel veel gecallt wordt (profiling maakt duidelijk dat ik >70% in deze functie vertoef) ben ik dus begonnen met het optimaliseren van de functie. Maar dat lukte ook niet echt, bitwise operators gaf geen snelheidsboost, en zelf OpenMP introduceren blijkt teveel overhead te geven (shared variable).

Momenteel ziet de code er al volgt uit:

C++:

// Compare two images
double compare(cairo_surface_t* inputSurface) const
{
        // Get the raw data of the given surface
        unsigned char* tempData1 = dataInputRGB24;
        unsigned char* tempData2 = cairo_image_surface_get_data(inputSurface);

        // Compare them
        int i, db, dg, dr;
        long int difference = 0;
        #pragma omp parallel for private(dr, dg, db) reduction(+:difference)
        for (i = 0; i < dataInputWidth*dataInputHeight*4; i+=4)
        {
                // RGBa
                db = tempData1[i] - tempData2[i];
                dg = tempData1[i+1] - tempData2[i+1];
                dr = tempData1[i+2] - tempData2[i+2];

                // Calculate difference
                difference += sqrt(dr*dr + dg*dg + db*db);
        }

        // Get resemblance
        double resemblance = dataInputWidth*dataInputHeight / double(difference);
        return resemblance;
}

Gprof:

code:

                0.00    1.15     204/204         EnvImage::fitness(DNA const&) const [4]
[5]     74.5    0.00    1.15     204         EnvImage::compare(_cairo_surface*) const [5]
                1.15    0.00     204/205         EnvImage::~EnvImage() [3]

Kan deze functie nog verder geoptimaliseerd worden? Of is het helemaal geen goede manier van afbeeldingen vergelijken, en bestaat er een beter/sneller alternatief?

Alvast bedankt
maleadt

[Maxxi]

Ik heb wel eens projecten gezien met beeldherkenning/vergelijking.

Meestal nemen ze daar meer dan 4 pixels als sample. Dat scheelt natuurlijk stuk.

Je forloop kan wel iets netter/sneller:


<code>
// Compare two images

int size = dataInputWidth*dataInputHeight*4;

for (i = 0; i < size ; i+=4)
{

</code>

Anders wordt elke keer opnieuw die berekening gedaan.

[Zoijar]

Beetje compiler ziet dat zelf wel hoor

Vergelijkings methode is een simpel MSE (mean-square-error) vergelijk. Vaak is dat prima, soms niet. Ligt eraan wat je wilt bereiken. Ik vergelijk zelf vaak in Lab space op delta E afstanden, dat is iets "perceptueler". Soms doen ze eerst een gaussian blur voor het vergelijk. Je kan ook naar dingen als VDP (visual differences predictor) or VEP (visual equivalence predictor) kijken, die zijn nog veel ingewikkelder met een model van hoe een mens iets ziet.

[Janoz]

Los van de gekozen vergelijkingsmethode denk ik dat de grootste kosten van je algoritme hem zitten in de wortel. Probeer je formules eens te herschrijven of je anders te formuleren zodat je niet voor elke pixel een wortel hoeft te doen.

Je zou bijvoorbeeld de wortel weg kunnen laten. Je score loopt dan wat sneller op wanneer de resultaten verder uit elkaar liggen, maar dat hoeft voor de scorings functie niet zo veel uit te maken.


Maar om even terug te komen op mijn eerste zin. Zou je wat meer kunnen vertellen over wat voor soort afbeeldingen je probeert te vergelijken? Cairo surfaces zeggen mij niet zo veel (tenzij je 1 of andere render lib bedoeld die ik met google tegen kwam, maar dan zegt dat nog weinig over wat voor plaatjes je nu daadwerkelijk vergelijkt). Je pixel per pixel vergelijking zorgt er nu namelijk wel voor dat een exat gelijk plaatje, dat een paar pixels verschoven is, absoluut niet als iets in de buurt gezien wordt.

[maleadt]

@Maxxi:
Bedankt voor die opmerking, maar zoals Zoijar zegt (en een snelle benchmark bevestigt) heeft de compiler dit door en cacht hij die waarde.

@Janoz:
Daar bedoelde ik inderdaad de Cairo rendering library mee. Maakt eigenlijk niet uit aangezien ze gewoon de raw pixel data aanbiedt, maar ik had het er maar bij gezet voor moest iemand weet hebben van een goed verborgen comparison functie (of iets dergelijks) in die library
En de wortel vermijden zorgt inderdaad voor een reductie van 1/3 in die routine, maar ik vrees dat de kwaliteit van m'n fitness functie er te zwaar gaat door achteruit gaan (en het is natuurlijk een nogal essentieel deel van een GA).

@Zoijar:
Delta E afstand lijkt me nog te doen, maar bij het opzoeken van die andere twee algoritmes begin ik te merken dat ik een complex (en dus vaak commerciëel) deel van de informatica betreed

EDIT; hmm, misschien iets leuk gevonden, eens kijken hoe het achter de schermen werkt.

[ Voor 6% gewijzigd door maleadt op 15-03-2009 19:41 ]

[TweakPino]

Je zou nog een for-lus om je for-lus kunnen knutselen, waarmee je OpenMP wel de gelegenheid geeft onafhankelijk te werken.

C++:

#define PARTS=2; 
long int difference_part[PARTS]; 
long int difference = 0;
int x;

#pragma omp parallel for 
for(x = 0; x < PARTS; x++)
{
    int i;
    difference_part[x] = 0;
    for (i = ((dataInputWidth*dataInputHeight*4)/PARTS) * x; i < ((dataInputWidth*dataInputHeight*4)/PARTS) * (x+1); i+=4)
    {
            // RGBa
            int db = tempData1[i] - tempData2[i];
            int dg = tempData1[i+1] - tempData2[i+1];
            int dr = tempData1[i+2] - tempData2[i+2];

            // Calculate difference
            difference_part[x] += sqrt(dr*dr + dg*dg + db*db);
    } 
}

for(x = 0; x < PARTS; x++)
{
    difference += difference_part[x];
}

[ Voor 7% gewijzigd door TweakPino op 15-03-2009 19:58 ]

[Soultaker]

Misschien kun je in plaats van OpenMP ook eens kijken naar de mogelijkheid om je algoritme met SSE instructies te implementeren. Je zou dan in ieder geval van vier 32-bit floating point getallen tegelijk de wortel moeten kunnen trekken. Misschien dat je de rest ook wat efficiënter kan uitvoeren.

[ Voor 12% gewijzigd door Soultaker op 15-03-2009 20:06 ]

[zoxzo]

Je zou afhankelijk van de GA generatie in eerste instantie minder pixels kunnen meenemen dwz iedere vierde bevoorbeeld en/of minder nauwkeurig dwz de wortel met slechts een paar termen van een taylorseries benadering berekenen.

Wat ook nog wel eens wil helpen is om de enkele loop in verschillende te splitsen:
1) eerst alle (gekwadrateerde) verschillen uit te rekenen, zonder de indices (i+3) over te slaan om de vectorisatie te helpen, daarna 2) de losse wortels, en als laatste 3) de sommatie, of een of andere variant hierop.

[maleadt]

@TweakPino:
Ha inderdaad, zo kan ik die reduction() laten vallen. Heb het nog geimplementeerd met omp_get_num_threads() ipv PARTS, en het geef enkele procentjes boost. De routine is echter nog altijd niet goed multithreaded, ik snap niet goed hoe dat komt. Misschien omdat ze eerder vaak gecallt wordt en niet zo processorintensief is.

@Soultaker:
SSE intrinsics zien er inderdaad heel interessant uit, maar ik kreeg ze ondanks een uurtje knutselen niet in gang. Bij gebrek aan low-level kennis is het allemaal net iets moeilijker, maar ik ga het zeker eens proberen. Als ik het goed begrijp kan ik enkel floats gebruiken als invoer voor SSE instructies?

@zoxzo:
Ik zou er de opbouw van mijn programma voor moeten veranderen, maar het is een interessant alternatief. Zeker als ik een complexere comparison functie ga gebruiken (want ik weet niet wat het beste werkt, een fitness() die snel werkt maar niet zo nauwkeurig is, of een hele trage fitness() die een heel nauwkeurig beeld van de huidige situatie geeft) is zo een nauwkeurigheidsevolutie best handig.

Dat is eigenlijk wel een belangrijk punt Is het beter om een fitness functie te hebben die misschien enkele honderden ms nodig heeft, maar daarvoor een heel nauwkeurige status kan neerzetten (hier: pixel-per-pixel analyse met een complex algorithme)? Of is het beter om een "goedkope" fitness-functie te gebruiken, in de hoop dat als een bepaalde fitness-waarde eens niet zo correct was, de runtime ervan zo kort was dat die mutatie wel gauw nog eens zal terugkomen (hier: maar 1 pixel per X evalueren en eenvoudig algo bezigen)?

[ Voor 0% gewijzigd door maleadt op 15-03-2009 23:02 . Reden: typo ]

[Matis]

Ik heb zelf ook een MSE functie geimplementeerd om ge-deinterlacte beelden te kunnen vergelijken met het orginele progressive beeld.

Om die MSE zo snel mogelijk te maken heb ik niet gebruikt voor sqrt(). Deze functie is zo traag en kost heel veel (onnodige) CPU tijd. Zoals hier boven ook al beschreven wordt.

Je zou ervoor kunnen kiezen om het hele process parralel te laten uitvoeren. Ik gebruik zelf CUDA (mocht je videokaart dat ondersteunen).

Mijn advies is om sws die sqrt() eruit te halen en om eens te kijken naar je variabele gebruik:

Immers gebruik je de int db, dg, dr voor het *bufferen* van de huidige pixel. Je zou dit ook direct kunnen doen; scheelt weer 3 keer schrijven naar een register.

[Nick The Heazk]

Janoz schreef op zondag 15 maart 2009 @ 19:18:
Los van de gekozen vergelijkingsmethode denk ik dat de grootste kosten van je algoritme hem zitten in de wortel. Probeer je formules eens te herschrijven of je anders te formuleren zodat je niet voor elke pixel een wortel hoeft te doen.

Deze opmerking is belangrijk. Gooi die wortel er maar uit. Met de som van de absolute waarden zul je waarschijnlijk een gelijkaardige performantie van je GA hebben (dat komt dan neer op een Manhattan distance).

Soultaker schreef op zondag 15 maart 2009 @ 20:06:
Misschien kun je in plaats van OpenMP ook eens kijken naar de mogelijkheid om je algoritme met SSE instructies te implementeren. Je zou dan in ieder geval van vier 32-bit floating point getallen tegelijk de wortel moeten kunnen trekken. Misschien dat je de rest ook wat efficiënter kan uitvoeren.

Ook dit is een belangrijke opmerking. Het verschil tussen de 3 kleurcomponenten kan prima met de SSE intrinsics van je compiler. Als je daar dan toch mee bezig bent, dan kan het interessant zijn om je lus wat te ontrollen, zodanig dat je minder jump instructies moet uitvoeren.

Verder kun je ook minder punten met elkaar vergelijken. Uiteindelijk kun je op die manier zelf bepalen hoeveel rekentijd je wilt besteden aan de compare functie. (Je geeft bijvoorbeeld een "skip factor" mee die aangeeft welke increment je gebruikt in je for lus).

[Soultaker]

maleadt schreef op zondag 15 maart 2009 @ 21:47:
@Soultaker:
SSE intrinsics zien er inderdaad heel interessant uit, maar ik kreeg ze ondanks een uurtje knutselen niet in gang. Bij gebrek aan low-level kennis is het allemaal net iets moeilijker, maar ik ga het zeker eens proberen. Als ik het goed begrijp kan ik enkel floats gebruiken als invoer voor SSE instructies?

Heel eenvoudig is het niet, maar je huidige optimalisaties zijn ook verre van eenvoudig, denk ik, dus wellicht is het een optie. Met MMX en SSE2 kun je trouwens wel rekenen met integers, maar uiteindelijk moet je toch naar float converteren om te kunnen worteltrekken.

Overigens lijkt zoxzo's suggestie om minder beeldpunten te vergelijken me ook wel zinnig, maar ik vraag me af of je dan niet gewoon je invoer terug kunt schalen (zodat al je plaatjes gewoon kleiner zijn) en dan kom je weer op hetzelfde algoritme uit. In principe zijn dit dus twee mogelijke verbeteringen die je los van elkaar kunt uitvoeren.

[Verwijderd]

Je gebruikt een 'unsigned char' om een kleurcomponent op te slaan. Dan is de grootste waarde die (dr*dr + dg*dg + db*db) kan hebben, 3 * 255^2 = 195020. De wortel hiervan is naar beneden afgerond 441, en past in een unsigned short. Dus kun je de sqrt() elimineren door gebruik te maken van een array 'unsigned short sqrt[195021]', waarbij je sqrt[i] compile-time gevuld hebt met de wortel uit i. In plaats van een dure call naar sqrt() doe je dan een aanzienlijk efficiëntere table-lookup.

Merk op dat je methode een resolutie van 0..441 heeft, en dat levert een tabel van 400K op - met mogelijk de nodige pageswapping. Je kunt overwegen genoegen te nemen met een wat lagere resolutie. Bijvoorbeeld: wanneer je (dr*dr + dg*dg + db*db) deelt door 4 (of 16), hetgeen voor int/unsigned een efficiënte operatie is, dan is je range 0..220 (of 0..110), en je tabel is nog maar 100K (of 25K). Wanneer deze tabel door de cpu gecached kan worden, levert dat een aanzienlijke winst op vergeleken met een grotere tabel.

Het downsamplen (terugschalen) van je images voor het vergelijken kan zinvol zijn. Er zijn grofweg twee methoden:

1. sampling van individuele pixels middels een selectie, bijvoorbeeld door het beschouwen van alleen even rijen en kolommen. Voordeel is dat dit efficiënt kan; nadeel is dat je driekwart van de informatie weglaat. Dat kan zeer goed werken, maar ook onaardige misinterpretaties opleveren bij bv zebrapad/schaakbord-achtige images;
2. sampling door het vervangen van groepen van pixels (bv blokken van 2x2) door een enkele pixel die de gemiddelde RGB waarde bevat. Voordeel tov. 1 is dat er minder informatie verloren gaat; nadeel is dat het rekenintensiever is. Indien je de images echter vaker dan een enkele keer vergelijkt, dan is deze manier te overwegen: je voert de downsampling dan maar een enkele keer uit, en bewaart het resultaat, dat je wanneer nodig hergebruikt.

Er zijn nog wel wat andere algorithmische optimalisaties denkbaar, maar die hangen sterk af van zaken als:

• is de kleur relevant of is vergelijking van de de vorm voldoende;
• wil je een 100% score halen, of is het voldoende om bv 80% van de gelijke images te vinden;
• hebben de 'beelden' van de images specifieke kenmerken waarop je eenvoudiger kunt filteren;
• hoevaak vergelijk je twee dezelfde images;
• hoeveel images vergelijk je: 10, 1K, 1M, ...

[maleadt]

@toaomatis:
sqrt heb ik ondertussen laten vallen, gaat de kwaliteit van mijn fitness() er niet teveel door achteruit

@Nick The Heazk & Soultaker:
tweede poging tot intrinsics is opnieuw op weinig uitgedraaid, maar ik begin de stof wel te snappen. Aangezien ik de sqrt heb laten vallen leken MMX intrinsics me beter, om zo minder registers te moeten gebruiken.

Aangezien mijn ingangsdata 1 byte is (unsigned char), en dus 4 bytes per pixel, leek me het eenvoudig om 1 __m64 per pixel te gebruiken (en zo wel 4 bytes te verkwisten, maar anders kan ik de substraction functies niet goed gebruiken).

Het begin van de "for" loop die alle pixels overloopt ziet er dan als volgt uit:

C++:

double compare_mmx(unsigned char* tempData1, unsigned char* tempData2)
{
    for (int i = 0; i < PIXELS*4; i+=4)
    {
        // Load both pixels
        __m64 mmxData1 = _mm_setr_pi16(tempData1[i], tempData1[i+1], tempData1[i+2], tempData1[i+3]); // 4x16bit {b, g, r, a}
        __m64 mmxData2 = _mm_setr_pi16(tempData2[i], tempData2[i+1], tempData2[i+2], tempData2[i+3]); // 4x16bit {b, g, r, a}
        
        // Substract
        __m64 mmxSubstraction = _mm_subs_pi16(mmxData1, mmxData2);  // 4x16bit {db, dg, dr, da}
        
        // Multiply
        __m64 mmxMult1 = _mm_mulhi_pi16(mmxSubstraction, mmxSubstraction);  // 2x32bit {db*db, dg*dg}
        __m64 mmxMult2 = _mm_mulhi_pi16(mmxSubstraction, mmxSubstraction);  // 2x32bit {dr*dr, da*da}
    }
    
    _mm_empty();
    return 0;
}

Klopt dit wat? Het probleem is echter de volgende stap, het optellen van de vier (of drie als ik "da" laat vallen) waarden in mmxMult1 en mmxMult2. De functie "_mm_add_pi32" toepassen op mmxMult1 en mmxMult2 levert me slechts een vector op met daarin "{db*db+dr*dr, dg*dg+da*da}", hoe kan ik die beide optellen? En hoe kan ik eventueel die "da" laten vallen, ik vind geen functies om een enkele bit uit een vector te halen?

@hij:
Leuke insteek, aan een statische tabel met voorberekende resultaten had ik eigenlijk echt niet gedacht. Zeker in combinatie met een lagere resolutie lookup is het een heel interessante aanpak
Dynamisch downsamplen zal ehcter enkel op de eenvoudige manier kunnen, gemiddelden van pixels berekenen wordt te duur aangezien bij elke vergelijking 1 van beide afbeeldingen volledig nieuw is en niet gecachet kan worden.

[zoxzo]

Het downsample idee was voornamelijk bedoelt om in de eerste iteraties te gebruiken, dwz wanneer je eigenlijk alleen geinteresseerd bent om snel uit je set van plaatjes een aantal relatief goede te kiezen. Dus moet je in latere iteraties ook steeds minder downsamplen en daarnaast kan je natuurlijk per iteratie voor het downsamplen alles een aantal pixels verschuiven zodat het verschil niet op dezelfde pixels blijft leunen.

@hij wijst er ook al op, maar als het even kan moet je kennis van hoe de afbeeldingen gemaakt worden of wat er in staat gebruiken om je vergelijkingsfunctie te optimaliseren.
Voorbeeldje:
Stel de doel-afbeelding bestaat uit een reeks overlappende rechthoeken van verschillend formaat en kleur etc en je wilt achterhalen hoe je dat met een minimaal aantal rechthoeken kan natekenen. Dus de GA genereert een aantal rechthoeken en bijbehorende kleuren. Vervolgens teken je
a) het plaatje en reken pixelwise het verschil uit.
Voordeel: werkt algemeen,
Nadeel: methode is relatief "dom" (maar desalniettemin verstandig om te testen).
b) van bovenop naar achter liggend en reken je alleen voor het zichtbare gedeelte van iedere rechthoek een verschil uit.
In het begin hoef je slechts voor de eerste rechthoek te optimaliseren, in een volgende GA iteratie hou je rekening met het daadwerkelijke kleurverschil, daarna stop je er een extra rechthoek achter en reken je dus alleen het verschil over de regio uit die daadwerkelijk verandert etc. etc.
Voordeel: per iteratie is het vergelijkingsgebied kleiner
Nadeel: allicht lastiger te implementeren, alhoewel je methode a) natuurlijk kan hergebruiken
--
Met de juiste byte gebaseerde constantes en bijpassende logical operator kan je natuurlijk het "a" of "da" gedeelte op nul zetten.
SSE3 heeft _mm_maddubs_pi16 die je op mmxSubtraction los kan laten om dan direct {db*db+dr*dr, dg*dg} te krijgen.

[ Voor 7% gewijzigd door zoxzo op 16-03-2009 22:36 . Reden: kleine toevoeging/verandering ]

[kunnen]

Probeer je code simpelweg eens met de intel compiler te compileren, en die te laten optimaliseren met SSE3, en alle standaard optimalisaties. Dit kan bij float-berekeningen tot wel 10x zo snel zijn als bijv. de microsoft compiler, in mijn ervaring.

[Nick The Heazk]

maleadt schreef op maandag 16 maart 2009 @ 21:15:
@Nick The Heazk & Soultaker:
tweede poging tot intrinsics is opnieuw op weinig uitgedraaid, maar ik begin de stof wel te snappen. Aangezien ik de sqrt heb laten vallen leken MMX intrinsics me beter, om zo minder registers te moeten gebruiken.

Aangezien mijn ingangsdata 1 byte is (unsigned char), en dus 4 bytes per pixel, leek me het eenvoudig om 1 __m64 per pixel te gebruiken (en zo wel 4 bytes te verkwisten, maar anders kan ik de substraction functies niet goed gebruiken).

Het begin van de "for" loop die alle pixels overloopt ziet er dan als volgt uit:
*knip code*

Klopt dit wat? Het probleem is echter de volgende stap, het optellen van de vier (of drie als ik "da" laat vallen) waarden in mmxMult1 en mmxMult2. De functie "_mm_add_pi32" toepassen op mmxMult1 en mmxMult2 levert me slechts een vector op met daarin "{db*db+dr*dr, dg*dg+da*da}", hoe kan ik die beide optellen? En hoe kan ik eventueel die "da" laten vallen, ik vind geen functies om een enkele bit uit een vector te halen?

Je bent op goede weg. Je kan die twee waarden optellen door de high en low 32-bits weg te schrijven naar twee gewone variabelen. Dan kun je de som op de normale wijze uitvoeren. In SSE3 heb je een HADD instructie waarmee je twee naast elkaar liggende words kunt optellen. Je hebt er tevens een instructie om de som te maken van alle words in je register en deze op te slaan in het eerste word. Waarschijnlijk heb je met MMX gelijkaardige instructies.

De da component kun je verwijderen door ze bij het inladen in je __m64 registers te vervangen door 0. Dus:

C++:

__m64 mmxData1 = _mm_setr_pi16(tempData1[i], tempData1[i+1], tempData1[i+2], 0); // 4x16bit {b, g, r, 0}

Let er overigens op dat de _mm_setr_pi16 jouw data in de omgekeerde volgorde in het register plaatst. Dit is mogelijk niet wat je wil. (Ik kan me vergissen, maar doorgaans staat die 'r' daar voor 'reverse').

[Matis]

Je heb t ook fucntie als memcmp waarme je hele geheugenblokken kunt vergelijken.

Ik weet niet precies wat er onder de motorkap gebeurt, maar je zou ervoor kunnen kiezen je data daar doorheen te kunnen gooien; Kijken wat voor variabelen je terugkrijgt (binnen het int domein).

Edit: heb geen ervaring met deze functie, maar ik weet van zijn bestaan af.

[ Voor 10% gewijzigd door Matis op 17-03-2009 09:14 . Reden: Extra tekst, Typo's ]

[maleadt]

@zoxzo:
Concreet voorbeeld, bedankt. Het valt inderdaad te zeggen dat mijn huidige manier van vergelijken niet al te slim is, maar ik betwijfel het of wat je voorstelt in mijn huidige opzet wel goed tot zijn recht komt. Aangezien het evolutionair gedeelte (mutate, crossover, ...) volledig losgekoppeld is van de omgeving (fitness, ...) kan ik maar moeilijk weten welke rechthoek er nu juist bijgekomen is. Of ik zou een kopie van het originele DNA moeten bijhouden en de verschillen berekenen, wat opnieuw neerkomt op de hele array overlopen en dus niet sneller zal zijn.
Maar tijdens dit schrijven komt er nog een alternatief in me op: de mutate() functie een soort delta parameter laten emitten... Je hebt iets aan het denken gezet, bedankt

@ThomasB:
memcpy was ook mijn eerste opzet, maar is niet nauwkeurig genoeg voor een evolutionair algoritme.

@Nick The Heazk:
En mijn code werkt min of meer toch

C++:

    // Total difference
    long int difference = 0;
    
    // Process all pixels
    for (int i = 0; i < PIXELS*4; i+=4)
    {
        // Load both pixels
        __m64 mmxPackedData1 = _mm_setr_pi16(tempData1[i], tempData1[i+1], tempData1[i+2], 0); // 4x16bit {b, g, r, 0}
        __m64 mmxPackedData2 = _mm_setr_pi16(tempData2[i], tempData2[i+1], tempData2[i+2], 0); // 4x16bit {b, g, r, 0}
        
        // Substract (packed arithmetics)
        __m64 mmxPackedSubstraction = _mm_subs_pi16(mmxPackedData1, mmxPackedData2);    // 4x16bit {db, dg, dr, 0}
        
        // Multiply and add (packed arithmetics)
        __m64 mmxPackedMAdd = _mm_madd_pi16(mmxPackedSubstraction, mmxPackedSubstraction);  // 2x32bit {db*db+dg*dg, dr*dr}
        
        // Add the result to the counter
        difference += _mm_extract_pi16(mmxPackedMAdd, 0) + _mm_extract_pi16(mmxPackedMAdd, 2);
    }

de HADD kan ik niet gebruiken, SSE3 lijkt hier niet te werken (Intel Linux x64).
Enkele vraagjes though
• • De code is niet echt sneller dan de normale methode. Ik denk dat het ligt aan het steeds converten van __m64 naar int. Als oplossing dacht ik om van "difference" een __m64 te maken, en de volledige 64bits (even groot als een long int) te gebruiken (niet packed dus). Maar hoe tel ik bij die variabele dan iets op? _mm_extract_* returnen allemaal ints, en gewone square brackets kan ik niet gebruiken om een __m64 met 32bit data uit een 64bit __mm64 te halen (difference += mmxPackedMAdd[0] + mmxPackedMAdd[1]). Of bestaan er intrinsics die met de volledige 64bits aan data werken (en waar ik dus overkijk)?
• _mm_madd_pi16 krijgt als input 2 4*16bit __m64 variabelen, en zou een 2*32 __m64 moeten teruggeven. Toch moet ik index 2 gebruiken om aan mijn data op de eerste plek te geraken? Schrijven brengt inzicht, natuurlijk moet dat aangezien _mm_extract_pi16 in stappen van 16 bit gaat. Maar daardoor verlies ik data, hoe extract ik 32 bit aan data uit deze vector (_mm_extract_pi32 bestaat niet in de MMX instructieset)?



Bedankt voor de reacties alvast!

[Soultaker]

Ik weet niet hoe goed je je code getest hebt, maar volgens mij gaat 't goed mis als het verschil tussen twee getallen groter dan 128 is; 128*128 is dan minstens 32768 en dan is je signed 16-bit integer negatief geworden. Als je bovendien twee van dat soort pixels optelt, past het geheel niet meer in een 16-bit word en mis je data.

Als je die laatste regel (optellen bij difference) zo herschrijft gaat 't wel goed:

C++:

        // Add the result to the counter
        difference += _mm_cvtsi64_si32(mmxPackedMAdd);
        difference += _mm_cvtsi64_si32(_mm_srli_si64(mmxPackedMAdd, 32));

...maar voor zover ik kon testen is de MMX implementatie zo niet sneller (eerder trager), dus het moet inderdaad wat anders.

maleadt schreef op dinsdag 17 maart 2009 @ 14:42:
• _mm_extract_* returnen allemaal ints, en gewone square brackets kan ik niet gebruiken om een __m64 met 32bit data uit een 64bit __mm64 te halen (difference += mmxPackedMAdd\[0] + mmxPackedMAdd\[1]). Of bestaan er intrinsics die met de volledige 64bits aan data werken (en waar ik dus overkijk)?
• [..] Maar daardoor verlies ik data, hoe extract ik 32 bit aan data uit deze vector (_mm_extract_pi32 bestaat niet in de MMX instructieset)?

Dit weet ik ook niet zo 1-2-3, maar het klopt wel dat niet alle denkbare SIMD instructies daadwerkelijk bestaan, dus je moet soms een beetje prutsen om efficiënt het gewenste resultaat te bereiken. (Zo is er wel een instructie om de laagste 32-bits te extracten, die ik in de code hierboven gebruiken, maar niet voor de hoogste.)

edit:
Is dat alpha-kanaal constant trouwens? Zo ja, dan kan het laden van pixels waarschijnlijk wat simpeler (en wat sneller).

[ Voor 3% gewijzigd door Soultaker op 18-03-2009 17:41 ]

[Soultaker]

Als ik het optellen ook in een MMX register doe en pixels in één keer inlaad (onder de aanname dat de alpha channel constant is, zodat ik die dus wel mee kan nemen) kom ik hier op uit:

C:

long long difference3(uint32_t *data1, uint32_t *data2)
{
    long long difference = 0; 
    int lo, hi;
    lo = 0;
    while (lo != PIXELS)
    {
        hi = lo + 16384;
        if (hi > PIXELS) hi = PIXELS;
        __m64 zero = _mm_setzero_si64();
        __m64 tmpsum = zero;
        do {
            __m64 x     = _mm_set_pi32(*data1++, *data2++);
            __m64 a     = _mm_unpacklo_pi8(x, zero);
            __m64 b     = _mm_unpackhi_pi8(x, zero);
            __m64 diff  = _mm_subs_pi16(a, b);
            __m64 sqsum = _mm_madd_pi16(diff, diff);
            tmpsum = _mm_add_pi32(tmpsum, sqsum);
        } while (++lo != hi);
        difference += (uint32_t)_mm_cvtsi64_si32(tmpsum);
        difference += (uint32_t)_mm_cvtsi64_si32(_mm_srli_si64(tmpsum, 32));
    }
    return difference;
}

Dit is toch wel significant (30% ofzo) sneller dan de puur-C++ implementatie.

[Exirion]

Mij is niet duidelijk wat je nou eigenlijk precies wil bereiken. Om wat voor afbeeldingen gaat het en wat doe je met de uitkomst van de 'gelijkeniscalculatie'? Je topic doet me eigenlijk denken aan dit topic:

\[C++] Objecten in foto herkennen

Als je op een abstracter niveau de gelijkenis tussen plaatjes wil bepalen dan moet je aan een heel ander soort algoritmen denken.

[Soultaker]

In aanvulling op de eerdere post: hetzelfde geintje kan ook met SSE2, dan kun je mooi vier paar pixels laden en twee paar pixels tegelijk verwerken:

C++:

long long difference4(__m128i *data1, __m128i *data2)
{
    assert(PIXELS%4 == 0);
    long long difference = 0;
    int lo = 0, hi = 0;
    while (lo != PIXELS)
    {
        hi += 16384;
        if (hi > PIXELS) hi = PIXELS;
        __m128i tmpsum = _mm_setzero_si128();
        do {
            __m128i x, y, a, b, diff, sqsum;

            x = _mm_load_si128(data1++);
            y = _mm_load_si128(data2++);

            a       = _mm_unpacklo_epi8(x, _mm_setzero_si128());
            b       = _mm_unpacklo_epi8(y, _mm_setzero_si128());
            diff    = _mm_subs_epi16(a, b);
            sqsum   = _mm_madd_epi16(diff, diff);
            tmpsum  = _mm_add_epi32(tmpsum, sqsum);

            a       = _mm_unpackhi_epi8(x, _mm_setzero_si128());
            b       = _mm_unpackhi_epi8(y, _mm_setzero_si128());
            diff    = _mm_subs_epi16(a, b);
            sqsum   = _mm_madd_epi16(diff, diff);
            tmpsum  = _mm_add_epi32(tmpsum, sqsum);
        } while ((lo += 4) != hi);
        difference += _mm_cvtsi128_si32(_mm_srli_si128(tmpsum,  0));
        difference += _mm_cvtsi128_si32(_mm_srli_si128(tmpsum,  4));
        difference += _mm_cvtsi128_si32(_mm_srli_si128(tmpsum,  8));
        difference += _mm_cvtsi128_si32(_mm_srli_si128(tmpsum, 12));
    }
    return difference;
}

Als je er nog meer uit wil halen kun je nog proberen om de kanalen onafhankelijk te verwerken (dan kun je makkelijker het alpha-kanaal overslaan) maar dan wordt de logica om data te laden weer wat ingewikkelder.

Eigenlijk is deze functie sowieso wel erg snel zo (aangezien je de wortel eruit gegooid hebt) dus waarschijnlijk is het zinniger je snelheidswinst elders te halen. De kans is groot dat de snelheid van je (cache) geheugen nu de grootste bottleneck is.

[Zoijar]

Beetje flauwe vraag tussendoor misschien, maar waarom moet het zo snel eigenlijk? Ik heb zelf ook wel eens zoiets gedaan met image vergelijkenm fouten te laten zien; daar moest ik zo ongeveer 32000 images vergelijken, en ik heb hem gewoon een keer 's avonds aangezet en 's ochtends was het klaar... stuk sneller dan dagen programmeren

[H!GHGuY]

Als je meer het menselijk oog wil benaderen voor de vergelijking dan kun je beter van RGB afstappen.
Yuv of een vergelijkbaar kleurensysteem is dan een stuk efficienter. Bovendien zal je waarschijnlijk meer belang hechten aan het verschil van het Y-kanaal dan aan die van de u/v kanalen.

[maleadt]

@Zoijar:
Ik schrijf dit programma gewoon omdat genetische algoritmes me interesseren, ik geef dus niet om een eenmalig resultaat maar probeer een mooi en efficiënt algoritme te hebben (waar ik echt geen dagen ga aan werken hoor). En daarbij is het wel leuk & leerrijk om een bepaalde functie eens ongeloofelijk hard optimaliseren.

@Exirion:
Genetisch algoritme dat een gegeven afbeelding benadert met een bepaald aantal polygonen. Heel m'n framework werkt, maar de evolutiesnelheid was me niet goed genoeg. In eerste instantie probeerde ik gewoon mijn fitness() functie (de zwaarste functie) te versnellen, maar o.a. uit dit topic leer ik dat het beter is een gesofisticeerder algoritme te pakken, maar dat dan niet domweg keer op keer de volledige afbeelding te laten verwerken maar intelligenter (bijvoorbeeld naargelang de generatie/mutatie slechts een bepaald deel heranalyseren).

@Soultaker 1:
Het ging inderdaad goed mis in eerste instantie, maar dat bedoelde ik met "min of meer" aangezien ik wist dat de extract_pi16 me dataverlies zou geven.

Interessante blokjes code ook. Maar (als ik het goed begrijp) kan ik dat 2e fragment toch niet zonder voorbewerking gebruiken, aangezien mijn data arrays unsigned char zijn (1 byte) en jij voor a en b 2 bytes uit de __m128 unpackt?
Het brengt mij wel op een idee, om via shifting operators een __m64 te misbruiken zodat ik mijn unsigned char[] direct aan een __m64 kan hangen. Daarbij kan ik de _mm_set_pi** vermijden, dewelke een dure operatie blijkt (om 4 bytes te setten genereert de compiler 40 instructies oid ).

[ Voor 0% gewijzigd door maleadt op 18-03-2009 20:50 . Reden: Eerst nog wat meer nadenken ]

[Soultaker]

maleadt schreef op woensdag 18 maart 2009 @ 20:46:
Interessante blokjes code ook. Maar (als ik het goed begrijp) kan ik dat 2e fragment toch niet zonder voorbewerking gebruiken, aangezien mijn data arrays unsigned char zijn (1 byte) en jij voor a en b 2 bytes uit de __m128 unpackt?

Juist wel (of ik begrijp je nu verkeerd); ik laad vier pixels op een rij in één XMM register (4x4 bytes = 16 bytes = 128 bits), en vervolgens unpack ik de onderste twee pixels van 8-bit bytes naar 16-bit words, voordat ik er verder mee ga rekenen. Daarna moet ik nog een keer hetzelfde doen voor de bovenste twee pixels.

Je moet het zo zien:

x =	A15	A14	A13	A12	A11	A10	A9	A8	A7	A6	A5	A4	A3	A2	A1	A0
y =	B15	B14	B13	B12	B11	B10	B9	B8	B7	B6	B5	B4	B3	B2	B1	B0
a =	0	A7	0	A6	0	A5	0	A4	0	A3	0	A2	0	A1	0	A0
b =	0	B7	0	B6	0	B5	0	B4	0	B3	0	B2	0	B1	0	B0
diff =	A7 - B7		A6 - B6		A5 - B5		A4 - B4		A3 - B3		A2 - B2		A1 - B1		A0 - B0
sqsum =	(A6 - B6)2 + (A7 - B7)2				(A4 - B4)2 + (A5 - B5)2				(A2 - B2)2 + (A3 - B3)2				(A0 - B0)2 + (A1 - B1)2

Hierna tel ik sqsum op bij tmpsum en herhaal ik de berekening voor de bovenste twee pixels van A en B (die staan nog steeds in x en y tenslotte).

tmpsum bevat de (tijdelijke) totaalsom en bestaat weer uit vier losse 32-bits dwords. Aangezien die kunnen overflowen, moet ik na zoveel iteraties even die waarden uit tmpsum halen en bij het totaal optellen. Dat gebeurt onderin de buitenste lus. De efficiëntie daarvan is niet interessant omdat die code niet vaak wordt uitgevoerd.

[voodooless]

Leuk allemaal , maar even een voorbeeldje om even in te haken op Exirion:

Neem een plaatje gevuld met een patroon van zwart en wit:

01
10

Alleen wat meer herhalingen. Neem nu een tweede plaatje om te vergelijken. Het tweede plaatje is ook een patroon:

01
10

't is een ander patroon, maar in essentie hetzelfde als het eerste, alleen een pixel verschoven in een willekeurige richting. Als je nu pixel per pixel de twee gaat vergelijken, dan zal de uitkomst zijn dat de plaatjes totaal niet op elkaar lijken. Ga je echter met het oog kijken, zul je vaststellen dat de plaatjes wel degelijk op elkaar lijken Nu is dus de vraag: wat wil je precies? Wat voor plaatjes ga je vergelijken, en op welke manier?

Dus voordat je dingen gaat versnellen moet je eerst maar eens gaan nadenken over hoe en wat lijkt me.

Anyway, het kan natuurlijk goed zijn dat de nu gebruikte methode voor de bedoelde toepassing prima werkt..