Discussion :

[sylvain.cool]

Bonjour,

J'ai une carte des altitudes de toute l'Europe que l'on peut apparenter en une image en niveau de gris. Celle-ci fait environ 5 milliards de pixels.
Mon but est de trouver les océans et les grands lac sur cette image. Mon approche est donc de trouver des régions suffisamment grande avec exactement la même altitude.
J'ai fait quelques essais et ça marche plutôt pas mal.

J'ai développé un petit algo tout simple moi-même. Je précise que je découpe d'abord l'image en zone plus petite pour ne pas avoir de buffer trop gros...

Code c++ :

Sélectionner tout - Visualiser dans une fenêtre à part

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
//Pour chaque pixel de la zone
for(int x=0; x<nSizeBuffX; x++)
{
	for(int y=0; y<nSizeBuffY; y++)
	{
		//On calcule la position du pixel dans le buffer
		int nPixPos = y*nSizeBuffX+x;
		//On init les paramètres
		int nRegionSize = 0;
		vRegion.clear();
		vPoints2Test.clear();
		//Si on n'est pas déjà passé par ce pixel, on regarde la taille de sa région
		if(!pAlreadyTestMap[nPixPos])
		{
			//On ajoute le premier point à la région
			POINT newPoint;
			newPoint.x = x;
			newPoint.y = y;
			vPoints2Test.push_back(newPoint);
			short int nCurAltitude = bufIn[nPixPos];
 
			//Tant qu'il reste des points à tester dans cette région
			while( vPoints2Test.size()>0 )
			{
				//On ajoute ce pixel à la région et on regarde si les pixels autour font partis de cette région
				nRegionSize += ComputeRegionSize(bufIn,nSizeBuffX,nSizeBuffY,nCurAltitude,																									pAlreadyTestMap,vRegion,vPoints2Test);
			}
 
			if( nRegionSize > fWaterMinSizePixels )
			{
				//Si la taille de la région est suffisament grande, on dit que c'est de l'eau
				SetWater(bufOut,vRegion,nCurAltitude,nSizeBuffX,nSizeBuffY);
			}
		}
	}
}

Code c++ :

Sélectionner tout - Visualiser dans une fenêtre à part

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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
int ComputeRegionSize(short int* bufIn, int &nSizeBuffX, int &nSizeBuffY, short int &nCurAltitude, bool* pAlreadyTestMap, vector<POINT> &vRegion, vector<POINT> &vPoints2Test)
{
	int x = vPoints2Test.begin()->x;
	int y = vPoints2Test.begin()->y;
	int nRegionSize = 1;
	int nPixPos = y*nSizeBuffX+x;
 
	//On ajoute le point courant à la region
	pAlreadyTestMap[nPixPos] = true;
	POINT newPoint;
	newPoint.x = x;
	newPoint.y = y;
	vRegion.push_back(newPoint);
 
	//On supprime le point courant des points à tester
	vPoints2Test.erase(vPoints2Test.begin());
 
	//On teste les points voisins
	int nNewX = x;
	int nNewY = y-1;
	if( nNewX>=0 && nNewX<nSizeBuffX && nNewY>=0 && nNewY<nSizeBuffY )
	{
		int nNewPixPos = nNewY*nSizeBuffX+nNewX;
		if( !pAlreadyTestMap[nNewPixPos] && nCurAltitude==bufIn[nNewPixPos] )
		{
			//Si jamais testé et même altitude, on les mets dans les points à tester
			newPoint.x = nNewX;
			newPoint.y = nNewY;
			vPoints2Test.push_back(newPoint);
			pAlreadyTestMap[nNewPixPos] = true;
		}
	}
 
	nNewX = x-1;
	nNewY = y;
	if( nNewX>=0 && nNewX<nSizeBuffX && nNewY>=0 && nNewY<nSizeBuffY )
	{
		int nNewPixPos = nNewY*nSizeBuffX+nNewX;
		if( !pAlreadyTestMap[nNewPixPos] && nCurAltitude==bufIn[nNewPixPos] )
		{
			//Si jamais testé et même altitude, on les mets dans les points à tester
			newPoint.x = nNewX;
			newPoint.y = nNewY;
			vPoints2Test.push_back(newPoint);
			pAlreadyTestMap[nNewPixPos] = true;
		}
	}
 
	nNewX = x;
	nNewY = y+1;
	if( nNewX>=0 && nNewX<nSizeBuffX && nNewY>=0 && nNewY<nSizeBuffY )
	{
		int nNewPixPos = nNewY*nSizeBuffX+nNewX;
		if( !pAlreadyTestMap[nNewPixPos] && nCurAltitude==bufIn[nNewPixPos] )
		{
			//Si jamais testé et même altitude, on les mets dans les points à tester
			newPoint.x = nNewX;
			newPoint.y = nNewY;
			vPoints2Test.push_back(newPoint);
			pAlreadyTestMap[nNewPixPos] = true;
		}
	}
 
	nNewX = x+1;
	nNewY = y;
	if( nNewX>=0 && nNewX<nSizeBuffX && nNewY>=0 && nNewY<nSizeBuffY )
	{
		int nNewPixPos = nNewY*nSizeBuffX+nNewX;
		if( !pAlreadyTestMap[nNewPixPos] && nCurAltitude==bufIn[nNewPixPos] )
		{
			//Si jamais testé et même altitude, on les mets dans les points à tester
			newPoint.x = nNewX;
			newPoint.y = nNewY;
			vPoints2Test.push_back(newPoint);
			pAlreadyTestMap[nNewPixPos] = true;
		}
	}
 
	return nRegionSize;
}

Maintenant, mon problème : il faut 2 jours à mon PC pour traiter tout ça !
Voyez vous quelque chose que je pourrais améliorer pour accélérer le traitement ? Le réduire à quelques heures serait pas mal...
J'aimerais rester sur mon approche, car je la trouve simple et surtout déjà implémentée :-), mais elle est surement assez limitée... donc si vous avez des idées pour tout changer et/ou utiliser des librairies comme opencv, pourquoi pas !

Merci d'avance,
Sylvain

[gpcbitnik38]

tu peux déjà vectoriser tes images tu va gagner pas mal en perf ^^
et je pense que tu peux aussi utiliser des image pyramidale:
tu crées une pyramide d'image (filtré avec un un filtre gaussien ) ça te permettra de trouver les grandes surfaces sur une petite image (la méditerranée la manche, l'atlantique etc... )
ensuite tu passe a l'image du niveau suivant ( plus résolu ) en prenant bien soin de supprimer les pixel déjà conservés.


après moi j'aurais tenté une autre approche :
utiliser un algo de clusterisation type kmeans ou meanshift, tu trouvera des implémentations optimiser et très rapide en ligne et avec un bon critère sur tes élévation ou ton gradient tu devrai pouvoir clusteriser sans trop de pb

[Smelk]

Tu peux peut être ne pas chercher à calculer la taille des régions dont l'altitude est trop élevée (si tu ne cherches que les océans) ?

[anapurna]

salut

si j'ai bien compris tu cherche les pixel de même couleur étant voisin entre eux
une fois trouvé ces voisins tu regarde si l'ensemble est assez étendu, si c'est effectivement le cas on considère que c'est une étendu d'eau

le plus simple aurait été d'utiliser l’algorithme de segmentation
te permettant de créer une liste de zone et après traiter les zones selon tes besoins

[sylvain.cool]

C'est bien ce que je fais Anapurna.

Une petite précision : ce que j'appelle mon "image" est une véritable carte d'élévation de l'Europe avec une précision de 1 ft en hauteur et de ~100m en horizontal. Chaque pixel a une élévation codée sur 2 octets.
Ceci fait que, hormis sur des zones d'océan ou lac, tout les pixels voisins auront une valeur différente.

Vu le nombre de pixel à traiter, je voyais un gros avantage à ma méthode, je ne passe qu'une seule fois sur chaque pixel.
Du coup, je ne vois pas bien ce que peuvent m'apporter des prétraitement alors que cet algorithme fonctionne déjà.

Je retiens l'idée du fonctionnement pyramidal car j'ai déjà la même image en définition réduite (en terme de nombre de pixel, pas de niveau d'élévation).
Au niveau des algos de segmentation, je vais regarder ce que propose gpcbitnik38 (kmeans ou meanshift). Mais est-ce que ce sera vraiment plus rapide... à voir !

Merci pour vos réponse en tout cas.

[anapurna]

salut

un liens qui t’intéresseras peut être ici

[Jipété]

Bonsoir,

Je suis peut-être hors-sujet mais, ne sachant pas d'où vient ta carte ni comment elle a été fabriquée, particulièrement concernant les espaces liquides, je me permets ce petit post :

[...] Mon approche est donc de trouver des régions suffisamment grandes avec exactement la même altitude.

[...] Une petite précision : ce que j'appelle mon "image" est une véritable carte d'élévation de l'Europe avec une précision de 1 ft en hauteur et de ~100m en horizontal.

[...] Ceci fait que, hormis sur des zones d'océan ou lac, tous les pixels voisins auront une valeur différente.

Dit autrement, à l'envers, ça donne Ceci fait que, juste sur des zones d'océan ou lac, tous les pixels voisins auront une valeur identique.

Et moi, là, je dis méfiance : as-tu déjà entendu parler des vagues scélérates ? Un truc de malade auquel personne ne croyait jusqu'à récemment, rare, totalement imprévisible et pouvant atteindre 30 m de haut !
Auquel cas tu seras loin, très loin, d'avoir des valeurs identiques, alors que c'est le même océan.

Regarde bien le dernier chapitre, "Observations & accidents notables", dans le lien que j'ai donné, il y en a encore eu une il y a tout juste un mois.

[pseudocode]

Maintenant, mon problème : il faut 2 jours à mon PC pour traiter tout ça !
Voyez vous quelque chose que je pourrais améliorer pour accélérer le traitement ? Le réduire à quelques heures serait pas mal...
J'aimerais rester sur mon approche, car je la trouve simple et surtout déjà implémentée :-), mais elle est surement assez limitée...

L'algo naif de "growing region" est assez lent par nature (exploration des pixels voisins, fifo, ...).

Regarde si un algo de "composantes connexes" ne serait pas plus rapide (par exemple l'algo "union-find". Il y a des exemple sur ce site).

[wiwaxia]

Bonjour,

Tu viens enfin de donner un facteur d'échelle

... Une petite précision : ce que j'appelle mon "image" est une véritable carte d'élévation de l'Europe avec une précision de 1 ft en hauteur et de ~100m en horizontal. Chaque pixel a une élévation codée sur 2 octets.
Ceci fait que, hormis sur des zones d'océan ou lac, tout les pixels voisins auront une valeur différente ...

Ce qui manque encore (ou qui m'a sans doute échappé), c'est la taille (ou l'aire) minimale des plans d'eau à mémoriser ... il serait intéressant de voir ce que donne la région de Scandinavie: la Finlande et la Russie voisine sont parsemées de milliers de lacs aux contours complexes, on peut comprendre que ton programme rame longtemps ...



... et rien n'a été dit au sujet des surfaces horizontales artificielles (stades, aéroports, parkings de supermarchés IKEA ...) Pour être juste, un écart minimal de 10 pixels correspond à une distance de l'ordre du kilomètre, et doit permettre d'éliminer les artefacts précédents.

Un test pas trop lourd pourrait concerner les 4 rangées de (n) pixels situées de part et d'autre du pixel central (x, y) et alignées avec lui, par le calcul des 4 sommes :
Sa = S_k=1^{à n}[Abs(P(x+k, y)- P(x,y))] , Sb = S_k=1^{à n}[Abs(P(x-k, y)- P(x,y))] ,
Sc = S_k=1^{à n}[Abs(P(x, y+k)- P(x,y))] , Sd = S_k=1^{à n}[Abs(P(x, y-k)- P(x,y))] ;
en comptant (1) pour chaque terme nul, tout point du plan d'eau serait caractérisé par une valeur comprise entre (1) et (4).

On pourrait y ajouter 4 termes diagonaux, afin de ne plus privilégier les directions parallèles aux axes:
Se = S_k=1^{à n}[Abs(P(x+k, y+k)- P(x,y))] , Sf = S_k=1^{à n}[Abs(P(x-k, y+k)- P(x,y))] ,
Sg = S_k=1^{à n}[Abs(P(x-k, y-k)- P(x,y))] , Sh = S_k=1^{à n}[Abs(P(x+k, y-k)- P(x,y))] ;
un point du plan d'eau présenterait alors une valeur comprise entre (1) et (8).

On respecterait ainsi les contours des rivages, en ne supprimant pas trop de détails.

@ Jipété: C'est juste, mais tu oublies l'aspect probabiliste: les rives du Léman n'ont pas la réputation d'être jonchées d'emballages de chocolat ou d'emmenthal, suite au déferlement de vagues scélérates . Si le repérage de ces dernières a été possible par satellite, c'est à la suite d'enquêtes méthodiques sur les zones d'accident (rives est de l'Afrique du Sud, notamment).
Ta remarque n'est cependant pas dépourvue de pertinence: en cas d'épisode météorologique agité, les écarts de pression atmosphérique peuvent produire des dénivellations mesurables sur de grandes distances, de l'ordre de quelques centaines de kilomètres; mais comme les photos satellites sont prises par beau temps, les conséquences sont probablement limitées ...

[wiwaxia]

Cela m'avait échappé, mais les marées apparaissent-elles sur l'image ? Elles viennent de l'Est, et la montée du niveau de la Mer du Nord précède celui de l'Atlantique avec un décalage certain.
D'où la nécessité d'un critère moins strict, tolérant sur de grandes distances des dénivellations liées aux amplitudes et déphasages des marées - consulter pour évaluation les tables appropriées; de plus, il ne faut pas oublier les erreurs accidentelles.

Admettre un dénivelé d'un pied (~ 30 cm) par kilomètre (soit 10 pixels) paraît raisonnable.

On reprendra donc le calcul des 4 premières sommes, pour (n) égal à 10, et en comptant (1) pour chaque terme au plus égal à 1 (ou 2).

Une sélection identique interviendra sur les quatre sommes diagonales:
Se = S_k=1^{à n₁}[Abs(P(x+k, y+k)- P(x,y))] , Sf = S_k=1^{à n₁}[Abs(P(x-k, y+k)- P(x,y))] ,
Sg = S_k=1^{à n₁}[Abs(P(x-k, y-k)- P(x,y))] , Sh = S_k=1^{à n₁}[Abs(P(x+k, y-k)- P(x,y))] ,
calculées cependant sur un plus petit nombre de termes:
n₁ = Round(n/2^1/2) = 7 afin de rester en-deçà de la même distance limite du pixel central (x, y).

PS [07/04/2017]: Puisqu'il suffit que l'un des huit termes précédents vérifie la condition: S_i < Limite pour que le point central appartienne à un plan d'eau, il est possible d'accélérer la recherche par l'arrêt des calculs successifs de (S₁, ... , S₈) dès que la condition en cause est satisfaite, au moyen d'une énumération conditionnelle du type: REPEAT ... UNTIL .

[sylvain.cool]

La finalité est d'afficher cette carte des élévations en mode topographie, mais en représentant les océans et "grands lacs" en bleu.
J'ai décidé unilatéralement qu'un grand lac faisait plus de 10 km² (ça a du bon d'être son propre chef :-)).
Ma carte est issue de donnée publique de l'UE (ici), donc je pense qu'ils ont réfléchi à tout ou partie des inconvénient donné plus haut.
Et enfin, ce n'est pas un problème si l'ensemble de l'océan atlantique ne fait pas une seule et même région. Il peut très bien être découpé en plusieurs petits ensembles.

Et je vais regarder plus en détails certaines des solutions proposées, mais en passant par des images pyramidales et en release plutôt qu'en debug (pas con hein ), je suis déjà passé à 7h (au lieu de 2 jours initialement)!

[sylvain.cool]

Et je viens de m'apercevoir que ce qui prend le plus de temps finalement, c'est la lecture écriture sur le disque...
Donc je pense que pour l'instant mon algo fait l'affaire... à vérifier quand je testerai ma carte finale !

[sylvain.cool]

Bonjour tout le monde,

Pour l'instant, mes essais donnent de bons résultats.
Reste le temps de lecture/écriture que je trouve un peu long... je vais poster un autre thread dans une autre section.
Et quand j'aurais un peu de temps, je regarderais les autres méthode proposée.

Merci beaucoup en tout cas pour vos réponses !