Optimisation algorithme d'érosion

Marsh Posté le 21-11-2008 à 15:14:42

Bonjour,

Je travaille sur une application qui a pour but de suivre des entités (des fourmis) sur des vidéo. Pour extraire les fourmis des images de la vidéo, je pratique une différence par rapport à une image de référence sur laquelle il n'y a aucune entité. J'obtiens ainsi une image qui ressemble à ça :

Comme on peut le voir l'image binaire est très bruité, pour la nettoyer je voudrai passer par une érosion. L'implémentation de l'algorithme en lui-même ne pose pas vraiment de problème. Le fait est surtout que l'exécution est plutôt lente... En effet, actuellement, je met environ 1s pour traiter une image de 720x576 pxl. Mon souci c'est qu'a raison de 25/s ça va prendre des heures a traiter la vidéo de 25min...

Donc je me demandais, si par hasard vous ne verriez pas une idée pour optimiser mon algo :

Code :

boost::shared_array<bool> FrameAnalyser::erosion(boost::shared_array<bool> image) {
boost::shared_array<bool> erodedImg(new bool[videoPlayer->getHeight()*videoPlayer->getWidth()]);
// Traitement des cas particuliers
for(int i = 0 ; i < videoPlayer->getWidth() ; ++i) {
erodedImg[i] = false;
erodedImg[videoPlayer->getHeight()*(videoPlayer->getWidth()-2)+i] = false;
}
for(int i = 0 ; i < videoPlayer->getHeight() ; ++i) {
erodedImg[i*videoPlayer->getWidth()] = false;
erodedImg[(i+1)*videoPlayer->getWidth()-1] = false;
}
// erosion de l'image
for(int j = 1 ; j < videoPlayer->getHeight() - 1 ; ++j) {
for(int i = 1 ; i < videoPlayer->getWidth() - 1 ; ++i) {
erodedImg[j*videoPlayer->getWidth()+i] = image[j*videoPlayer->getWidth()+i-1]
&& image[j*videoPlayer->getWidth()+i]
&& image[j*videoPlayer->getWidth()+i+1]
&& image[(j-1)*videoPlayer->getWidth()+i]
&& image[(j+1)*videoPlayer->getWidth()+i];
}
}
return erodedImg;
}

Car là, à par aller voir les optimisation du côté de l'accès au cache... je vois pas trop ce que je peux faire...
Après je me dit que ce sera pas l'algorithme le plus couteux, donc si il y a pas d'autres solutions que de traiter l'accès au cache. Je ferai mieux de laisser tomber et d'aller voir la suite (je perdrai moins de mon temps...)

Merci

[edit] si vous avez d'autres idées qu'une érosion pour nettoyer l'image je suis preneur aussi

Message édité par Amonchakai le 21-11-2008 à 15:20:31

Reply

Marsh Posté le 21-11-2008 à 15:14:42

Reply

Marsh Posté le 21-11-2008 à 15:41:39

Déjà, tu pourrais peut-être multithreader le traitement non !? Chaque thread traite une image, et donc si multi-core ou autre, le traitement dans son ensemble peut-être largement accéléré !?

---------------
Kao ..98 - Uplay (R6S) : kao98.7.62x39 - Origin (BF4, BF1) : kntkao98

Reply

Marsh Posté le 21-11-2008 à 16:19:27

stop, stop Stop.
Avant d'aller plus loin, il faut optimiser ton algo pr etre cache friendly, ce qui represente des gains de l'ordre de 5à 10.

Donc, faut maximiser la localité spatiale en découpant ton image en bloc carré de taille approximativement égale à celle d'un bloc de cache line (genre 32x32 ou 16x16) et appliquer ton algo dessus.

En gros tu te retrouve avec 4 nid de boucle : un sur les tuile en x, en y et pour chaque tuile en x et y.

Aprés, l'accès des tableaux 2D est mieux si ton tableau est 2D au sens des NRC, ce que ne supporte pas forcement shared_array. J'ai deja poster 5 fois le code, je te laisse faire la recherche.

Ensuite, vectorise en utilisant SSE2

La déjà tu devrais avoir gagner en 10 et 40.

Si ca va pas assez vite , alors oui : boost::thread.

PS :aussi, on evite d'allouer de la mémoire dans les fonctions de calcul On passe la sortie en paramètre references
PS2 : factorise tes calculs d'adresse aussi

Message cité 1 fois
Message édité par Joel F le 21-11-2008 à 16:23:28

Reply

Marsh Posté le 21-11-2008 à 21:19:23

qqun a dit boost::gil ??

Reply

Marsh Posté le 21-11-2008 à 21:33:57

non car les perfs sont atroces.

Message cité 1 fois

Reply

Marsh Posté le 21-11-2008 à 22:40:50

Okay, super, je vais aller voir tout ça.

Merci

Message édité par Amonchakai le 21-11-2008 à 22:41:05

Reply

Marsh Posté le 22-11-2008 à 02:34:49

Joel F a écrit :

stop, stop Stop.
Avant d'aller plus loin, il faut optimiser ton algo pr etre cache friendly, ce qui represente des gains de l'ordre de 5à 10.

Donc, faut maximiser la localité spatiale en découpant ton image en bloc carré de taille approximativement égale à celle d'un bloc de cache line (genre 32x32 ou 16x16) et appliquer ton algo dessus.

En gros tu te retrouve avec 4 nid de boucle : un sur les tuile en x, en y et pour chaque tuile en x et y.

Aprés, l'accès des tableaux 2D est mieux si ton tableau est 2D au sens des NRC, ce que ne supporte pas forcement shared_array. J'ai deja poster 5 fois le code, je te laisse faire la recherche.

Ensuite, vectorise en utilisant SSE2

La déjà tu devrais avoir gagner en 10 et 40.

Si ca va pas assez vite , alors oui : boost::thread.

PS :aussi, on evite d'allouer de la mémoire dans les fonctions de calcul On passe la sortie en paramètre references
PS2 : factorise tes calculs d'adresse aussi

Avant d'aller aussi là, on en revient à la définition des fonctions Height et Width, est-ce inliné ou pas, et surtout est-ce que ça va être relu à chaque appel ou pas.
Et surtout la traditionnelle histoire sur la boucleavec gcc, est-ce qu'il va optimisé le j*h ou bien est-ce qu'il va falloir préférer du j += h; ...

ma question bonus ça serait: est-ce que si on passe sur du stockage bit à bit on peut y gagner en localité (et bien déroulant les tests/sets des 8 bits) ou bien est-ce que les opérations sur les bits sont plus lentes et mieux vaut lire un bon vieil octet ? t'as déjà comparé ?

Message cité 1 fois

Reply

Marsh Posté le 22-11-2008 à 11:09:43

Joel F a écrit :

non car les perfs sont atroces.

Si t as un lien vers un article/bench/post qui parle de ca je suis interesse.

Reply

Marsh Posté le 22-11-2008 à 11:28:56

la mailing de boost contient pas de sujet à ce propos.
gil c'est bien pr les io. le reste ...

aprés moi, vu que c'est mon boulot d'implanter ce genre de truc, j'ai regardé le source de gil. C'est pas cahce-friendly, ca itére avec de siterateurs au lieu d'index ... bref.

Reply

Marsh Posté le 22-11-2008 à 11:29:41

Taz a écrit :

Avant d'aller aussi là, on en revient à la définition des fonctions Height et Width, est-ce inliné ou pas, et surtout est-ce que ça va être relu à chaque appel ou pas.

oui d'ou mon PS

Taz a écrit :

Et surtout la traditionnelle histoire sur la boucleavec gcc, est-ce qu'il va optimisé le j*h ou bien est-ce qu'il va falloir préférer du j += h; ...

Plus trop maintenant. avec gcc 4.3+, la différence est de l'ordre du bruit des mesures sur Intel.

Taz a écrit :

ma question bonus ça serait: est-ce que si on passe sur du stockage bit à bit on peut y gagner en localité (et bien déroulant les tests/sets des 8 bits) ou bien est-ce que les opérations sur les bits sont plus lentes et mieux vaut lire un bon vieil octet ? t'as déjà comparé ?

ca sert à rien, parole d'expert . Je serais aussi enclin à stocker des uchar plutot que des bool dont la taille varie en fonction du compilo.

Message édité par Joel F le 22-11-2008 à 11:35:08

Reply

Marsh Posté le 22-11-2008 à 11:29:41

Reply

Marsh Posté le 22-11-2008 à 11:31:09

Je viens de tester un utilisant les trucs que je file en TP. J'arrive à 13 cycles par point pr une erosion sans optimisation. Soit, 2.8ms pour une image 720*576
sur un PC à 2GHz. Soit un x480 de mieux que le code original. Je subodore que les fonctions d'acces au taille sont pas inline et que tu meurt en allouant ton image résultat dans ta fonction.

Le code :

Code :

// Mesure cycles pour Intel & gcc
static inline volatile double read_cycles()
{
uint32_t hi, lo;
__asm __volatile ("rdtsc" : "=a" (lo), "=d" (hi));
return double((long long)hi << 32 | lo);
}
// Gestion mémoire
template<class T> static inline T** alloc( size_t h, size_t w )
{
T **m;
m = (T**)malloc( h*sizeof(T*) );
m[0] = (T*)malloc( h*w*sizeof(T) );
for(size_t i=1;i<h;i++) m[i]=m[i-1]+w;
return m;
}
template<class T> static inline void release( T** ptr )
{
if(ptr) free(ptr[0]);
if(ptr) free(ptr);
}
// Calcul
void erosion(unsigned char** I,unsigned char** O,int i0,int i1,int j0,int j1)
{
unsigned char r,u,ml,mm,mr,d;
for(int i=i0+1;i<=i1-1;++i)
for(int j=j0+1;j<=j1-1;++j)
{
// Element structurant scalarisé
u = I[i-1][j];
ml = I[i][j-1]; mm = I[i][j]; mr = I[i][j+1];
d = I[i+1][j];
// Evaluation
r = (u && d) && (ml && mr && mm);
// Stockage
O[i][j] = r;
}
}
// Appel
int main()
{
static const size_t h = 720, w = 576, iter = 1000;
unsigned char** input = alloc<unsigned char>(h,w);
unsigned char** output = alloc<unsigned char>(h,w);
double cs,c0;
c0 = read_cycles();
for(size_t i=0;i<iter;++i) erosion(input,output,0,h-1,0,w-1);
cs = read_cycles() - c0;
cout << "cpp scalaires : " << cs/(iter*h*w) << endl;
cout << "time scalar : " << cs/(iter*2e9)*1000 << "ms" << endl;
release(input);
release(output);
return 0;
}

Citation :

Execution :
cpp scalaires : 13.1564
time scalar : 2.72811ms

Message édité par Joel F le 22-11-2008 à 14:29:16

Reply

Marsh Posté le 22-11-2008 à 16:47:00

Merci pour l'info.

Un petit site pas mal en passant:
http://www.cellperformance.com/mike_acton/

avec certains articles chouettes comme http://www.cellperformance.com/mik [...] _keyw.html

et le dernier papier de Drepper (GNU libc entre autres) avec un bon chapitre 6 sur le cache avec un petit exemple
http://people.redhat.com/drepper/cpumemory.pdf

Message édité par Taz le 22-11-2008 à 16:51:37

Reply

Marsh Posté le 22-11-2008 à 17:15:36

mike acton est pas mal.
par contre, en general, le cache suffit pr ce genre de truc. les optim de fsb ou de vram : bof ^^

Deja si tout le mode faisait du LCS patten et du blocking dans ces nids d eboucles :E ca serait pas mal

Message édité par Joel F le 22-11-2008 à 17:19:50

Reply

Marsh Posté le 22-11-2008 à 20:37:16

ça reste un bon article pour avoir une vision d'ensemble

Reply

Marsh Posté le 22-11-2008 à 21:05:28

C'ets una rticle du genre que j'utilise pr mes cours d'archi. Va fallori que je me mettes à jour

Reply

Marsh Posté le 29-06-2009 à 12:56:41

Bonjour,
Je ressors ce thread car j'ai eu récemment à réimplémenter cet algo et j'avoue être assez perplexe sur les gains que j'obtiens après une optimisation vis a vis du cache...

voici mon algo:

Operator.hpp

Code :

static const size_t BLOCK = 32;
template<class T> static inline T** alloc( size_t h, size_t w ) {
T **m;
m = (T**)malloc( h*sizeof(T*) );
m[0] = (T*)malloc( h*w*sizeof(T) );
for(size_t i=1;i<h;i++) m[i]=m[i-1]+w;
return m;
}
template<class T> static inline void release( T** ptr ) {
if(ptr) free(ptr[0]);
if(ptr) free(ptr);
}
// ----------------------------------------------------------------------------------------------------------------------
// apply a block based operator to a picture
template<typename Operator, int radius> static inline
void apply(const unsigned char * const* in, unsigned char **out, int width, int height) {
size_t lastVBlock = height / BLOCK;
size_t lastHBlock = width / BLOCK;
size_t posPrevLastVBlock = (lastVBlock-1)*BLOCK;
size_t posPrevLastHBlock = (lastHBlock-1)*BLOCK;
BOOST_STATIC_ASSERT(radius < BLOCK);
for(size_t j = 1 ; j < lastVBlock - 1 ; ++j) {
for(size_t i = 1 ; i < lastHBlock - 1 ; ++i) {
size_t iBlock = i*BLOCK;
size_t jBlock = j*BLOCK;
Operator::template apply<radius>(in, out, iBlock, jBlock, iBlock+BLOCK, jBlock+BLOCK);
}
}
for(size_t j = 1 ; j < lastVBlock - 1 ; ++j) {
size_t jBlock = j*BLOCK;
Operator::template apply<radius>(in, out, radius, jBlock, BLOCK, jBlock+BLOCK);
Operator::template apply<radius>(in, out, posPrevLastHBlock, jBlock, width-radius, jBlock+BLOCK);
}
for(size_t i = 1 ; i < lastHBlock - 1 ; ++i) {
size_t iBlock = i*BLOCK;
Operator::template apply<radius>(in, out, iBlock, radius, iBlock+BLOCK, BLOCK);
Operator::template apply<radius>(in, out, iBlock, posPrevLastVBlock, iBlock+BLOCK, height-radius);
}
Operator::template apply<radius>(in, out, radius, radius, BLOCK, BLOCK);
Operator::template apply<radius>(in, out, radius, posPrevLastVBlock, BLOCK, height-radius);
Operator::template apply<radius>(in, out, posPrevLastHBlock, posPrevLastVBlock, width-radius, height-radius);
Operator::template apply<radius>(in, out, posPrevLastHBlock, radius, width-radius, BLOCK);
// Operator::template apply<radius>(in, out, radius, radius, width-radius, height-radius);
}
// ----------------------------------------------------------------------------------------------------------------------
// Erosion
struct ErodeByBlock {
template<int size> static inline
void apply(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height);
} ;
template<int radius>
inline void ErodeByBlock::apply(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height) {
for(int j = yOff ; j < height ; ++j) {
for(int i = xOff ; i < width ; ++i) {
out[j][i] = 255;
for(int jj = -radius ; jj < radius ; ++jj) {
for(int ii = -radius ; ii < radius ; ++ii) {
if(std::abs(ii)+std::abs(jj) <= radius)
out[j][i] = std::min(out[j][i], in[j+jj][i+ii]);
}
}
}
}
}
template<>
inline void ErodeByBlock::apply<1>(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height) {
for(int j = yOff ; j < height ; ++j) {
for(int i = xOff ; i < width ; ++i) {
unsigned char xnm0 = std::min(in[j-1][i], in[j][i-1]);
unsigned char xnm1 = std::min(in[j][i], in[j+1][i]);
xnm0 = std::min(xnm0, in[j][i+1]);
out[j][i] = std::min(xnm0, xnm1);
}
}
}
template<int radius> static inline
void erosion(const unsigned char * const* in, unsigned char **out, int width, int height) {
apply< ErodeByBlock, radius >(in, out, width, height);
}
// ----------------------------------------------------------------------------------------------------------------------
// Dilatation
struct DilateByblock;
template<int radius> static inline
void dilatation(const unsigned char * const* in, unsigned char **out, int width, int height) {
apply< DilateByblock, radius >(in, out, width, height);
}
struct DilateByblock {
template<int size> static inline
void apply(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height);
} ;
template<int radius>
inline void DilateByblock::apply(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height) {
for(int j = yOff ; j < height ; ++j) {
for(int i = xOff ; i < width ; ++i) {
out[j][i] = 0;
for(int jj = -radius ; jj < radius ; ++jj) {
for(int ii = -radius ; ii < radius ; ++ii) {
if(std::abs(ii)+std::abs(jj) < radius)
out[j][i] = std::max(out[j][i], in[j+jj][i+ii]);
}
}
}
}
}
template<>
inline void DilateByblock::apply<1>(const unsigned char * const* in, unsigned char **out, int xOff, int yOff, int width, int height) {
for(int j = yOff ; j < height ; ++j) {
for(int i = xOff ; i < width ; ++i) {
unsigned char xnm0 = std::max(in[j-1][i], in[j][i-1]);
unsigned char xnm1 = std::max(in[j][i], in[j+1][i]);
xnm0 = std::max(xnm0, in[j][i+1]);
out[j][i] = std::max(xnm0, xnm1);
}
}
}

main.cpp

Code :

static inline volatile double read_cycles() {
uint32_t hi, lo;
__asm __volatile ("rdtsc" : "=a" (lo), "=d" (hi));
return double((long long)hi << 32 | lo);
}
int main() {
size_t width = 736;
size_t height = 595;
size_t iteration = 1000;
double cs,c0;
c0 = read_cycles();
unsigned char **image = alloc<unsigned char>(width, height);
unsigned char **out = alloc<unsigned char>(width, height);
for(size_t i = 0 ; i < iteration ; ++i) {
erosion<1>(image, out, width, height);
dilatation<1>(image, out, width, height);
}
cs = read_cycles() - c0;
std::cout << "cpp scalaires : " << cs/(iteration*height*width) << std::endl;
std::cout << "time scalar : " << cs/(iteration*2e9)*1000 << "ms" << std::endl;
release<unsigned char>(image);
release<unsigned char>(out);
return 0;
}

Bilan, la version cache friendly ne m'apporte quasi aucune amélioration sur les vidéos sur lesquelles je travaille (en CIF)

Citation :

défaut:
cpp scalaires : 34.6357
time scalar : 1.75561ms

Optimisé:
cpp scalaires : 34.6327
time scalar : 1.75546ms

Si j'augmente la taille (du 4CIF comme dans l'exemple posté)

Citation :

défaut:
cpp scalaires : 76.1364
time scalar : 16.6708ms

optimisé:
cpp scalaires : 63.462
time scalar : 13.8956ms

C'est toujours pas folichon folichon comme amélioration... Bon, niveau perf, honnêtement, ça me va... je me demandais juste si vous aviez une explication du si peu de gain que j'obtiens... Ou, peut-être ai-je encore fait une grosse erreur...

Merci

[edit] qq infos sur ma config, je suis sous Linux avec gcc 4.4.0-3, le tout compilé en O2

Message édité par Amonchakai le 29-06-2009 à 13:06:42

Reply

Marsh Posté le 29-06-2009 à 13:11:19

euh y a qd meme 13 CYCLES de gagné par points
Sinon -O3 aussi car pas sur que -O2 inline à fond.
Tu gagneras aussi à pas faire de if(std::abs()) dans tes boucles, t'as pas beosin de ça (cf vieux post)

Pour aller plus loin, SIMDisation ou bien allocation alignée sur la cache-line

Message édité par Joel F le 29-06-2009 à 13:13:52

Reply

Marsh Posté le 29-06-2009 à 13:23:10

Oui, j'ai essayé avec SSE2.. mais très mal fait

Ce qu'il y a c'est qu'avec mon érosion en 4 connexité avec un rayon de 1, ça fait pas beaucoup nombres donc mes vecteurs sont pas suffisamment remplit...
Du coup, j'avais cassé complètement les perfs (le compilo & le proc se débrouillent mieux sans moi)

Je pense qu'il faudrait que je déploie les boucles pour remplir mieux les vecteurs...

Reply

Marsh Posté le 29-06-2009 à 13:45:30

Je pense surtout que ta convolution etait mal vectorisé
Elimine le if deja tu verras, ca ira mieux.

Pour ref:
http://www.hipeac.net/system/files [...] evised.pdf

Saute le passage sur le CELL, les optimsiations présentées sont valides sur x86 et cie.

Message édité par Joel F le 29-06-2009 à 13:45:50

Reply

Marsh Posté le 29-06-2009 à 13:54:01

Le if, c'est bon je ne l'ai pas... Il est juste présent dans la version générique, j'ai fait une spécialisation pour un rayon de 1 (et 2 et 3)

Et merci, pour le document je vais étudier ça

[edit]
Après dans les versions générique je pense que j'ai vraiment besoin de faire les if(abs()), car je souhaite assurer une distance en 4 connexité

Message édité par Amonchakai le 29-06-2009 à 14:00:25

Reply

Marsh Posté le 29-06-2009 à 13:55:46

le if() est inutile tt le temps

Reply

Marsh Posté le 29-06-2009 à 14:01:15

ha... je vais voir

[edit]
On est d'accord, je sors pas de la zone allouée. Mais par contre, je passe en 8 connexité si j'enlève le if non ?

Message cité 1 fois
Message édité par Amonchakai le 29-06-2009 à 14:08:35

Reply

Marsh Posté le 29-06-2009 à 16:15:49

Amonchakai a écrit :

ha... je vais voir

[edit]
On est d'accord, je sors pas de la zone allouée. Mais par contre, je passe en 8 connexité si j'enlève le if non ?

Oui certes.

Reply

Marsh Posté le 03-07-2009 à 15:53:04

Hello,

Bon, Il s'agit juste un post un peu inutile...
Mais je voulais juste redire merci au sujet de l'article qui a été mis en lien précédemment. Il y a plein de choses intéressantes dedans (il s'avère que j'ai beaucoup de boulot sur mes convolutions par des filtres 2D). Mais, il contient plein de bon principes (notament des réponses a des questions que je me posais sur l'apport du regroupement des calculs...)

Donc merci encore

# Avec SSE2, j'arrive maintenant a gratter encore quelques 40%. donc c'est pas mal du tout (et ça a eu le mérite d'être instructif)

Message édité par Amonchakai le 03-07-2009 à 16:01:59

Reply

Marsh Posté le 27-07-2009 à 18:16:27

Amonchakai, plutot que de te lancer dans de l'optimisation bas-niveau direct, re-reflechis peut-etre au probleme initial.
De maniere generale, l'erosion/dilatation, c'est du rafistolage, ça se generalise tres mal et c'est bourrain.

T'as regarde du cote du water-shedding ?

Si tu tiens réellement a jouer avec dilatations/erosions, lance une convolution avec un filtre gaussien a la con, tu reduira le nombre d'iterations necessaires pour bien decouper tes fourmis :sweat:

Reply

Marsh Posté le 28-07-2009 à 11:17:29

water-shedding c'est quand meme la massue bien lente pour ce genre de probleme.
Par contre oui, gaussian blur avant.

Message édité par Joel F le 28-07-2009 à 11:17:41

Reply

Marsh Posté le 28-07-2009 à 11:37:22

Hello,

Merci pour cette remarque. Mais en réalité j'avais réglé le problème en générant une image de référence sans fourmi et soustrait cette image aux images contenant les fourmis + un seuillage.
Ca marche plutôt bien.

Pour ce qui est de l'érosion/dilatation j'en ais juste eu besoin récemment pour un algo d'ouverture/fermeture par reconstruction... Et j'ai lu des papiers qui décrivaient une implémentation efficace à base d'érosion et de dilatation...

enfin, merci pour l'aide

Reply

Marsh Posté le

Reply

Optimisation algorithme d'érosion

Sujets relatifs:

Leave a Replay