Resultados | Exercícios 6 | Alexandre Felipe (oalexandrefelipe)

activity_oalexandrefelipe_6_twsm - Tie-zone Watershed from labeled markers

1. Definição da função a ser feita e testada

Implementar a transformada tie-zone watershed por marcadores rotulados similar à ia870:iatz, porém com a opção onde os marcadores já entram rotulados e o resultado é uma imagem rotulada e não as linhas do watershed.

O objetivo desta implementação é dar uma experiência sobre a implementação do watershed. A parte fundamental é a fila hierárquica e o desempenho do algoritmo é diretamente relacionada à eficiência desta fila hierárquica.

Sinta-se à vontade para fazer as otimizações que achar conveniente.

 1 import numpy as np
 2 import ia870;
 3 import heapq; # Heap queue built in to python 2.3
 4 
 5 def dijkstra(f, m, Bc = None):
 6   # Index of the items to be visited
 7   Bci = zip(*np.nonzero(Bc));
 8   # Auxiliary Variables
 9   lbl = m.copy(); # Labels
10   f_inf = np.max(f) + 1;
11   costs = np.tile(f_inf, f.shape);
12   cost2 = np.zeros(f.shape);
13   costs[m != 0] = f[m != 0]; # Initial costs
14   done = np.zeros(f.shape, np.bool);
15   # Insert initial items in a heap queue
16   mi, mj = m.nonzero();
17   mi, mj = mi.ravel(), mj.ravel();
18   hq0 = zip(f[mi, mj], np.zeros(mi.shape), mi, mj);
19   heapq.heapify(hq0);
20   pos = len(mi);
21   while(len(hq0) != 0):
22     # Source cost and coordinates.
23     _cost, _cost2, _i, _j = heapq.heappop(hq0);
24     done[_i, _j] = True;
25     # The cost used in the queue is the same
26     # cost saved in the cost array.
27     assert(_cost == costs[_i, _j]);
28     for i,j in Bci:
29       # destination pixel
30       di, dj = _i + i - 1, j + _j - 1;
31       if(np.all((di >= 0, di < f.shape[0],
32                  dj >= 0, dj < f.shape[1]))):
33         if done[di, dj]:
34           continue; # This is not necessary
35         # Destination pixel is inside the image.
36         f_dest = f[di, dj];
37         n_cost = max(_cost, f_dest);
38         if(n_cost == costs[di, dj]):
39           if(lbl[di,dj] != lbl[_i,_j]):
40             # lbl[di, dj] = 0; # What I initially implemented
41             # This whould increase the tie zone region.
42 
43             if(n_cost == costs[_i, _j]):
44               # If it is a plateau, uses the distance
45               # to the border of the plateau to decide
46               # whether it is a tie or not.
47               if(cost2[di, dj] == cost2[_i, _j] + 1):
48                 lbl[di,dj] = 0; # Tie zone
49             else:
50               # It is not in a plateau, the cost is increasing
51               # In this case we know that the path used to
52               # reach the point (di, dj) ended with a plateau
53               # or is the only pixel having a value of dissimilarity
54               # so large.
55               lbl[di, dj] = 0; # Tie zone.
56         elif(n_cost < costs[di, dj]):
57           # We have just reached to the destination point
58           # by a path whith lower cost than the previously
59           # computed value.
60           # Due to the heap queue implementation, this code
61           # will run only once each pixel.
62 
63           # Make sure each point is visited at most once.
64           assert(costs[di,dj] == f_inf);
65           lbl[di, dj] = lbl[_i, _j];
66           costs[di, dj] = n_cost;
67           if(n_cost == costs[_i, _j]):
68             # Distance to the border of the plateau.
69             cost2[di, dj] = cost2[_i, _j] + 1;
70           heapq.heappush(hq0, (n_cost, cost2[di, dj], di, dj));
71 
72   assert(np.all(done)); # Check whether all pixel was calculated.
73   return lbl, costs, cost2;
74 
75 
76 def twsm(f,m,Bc):
77   lbl, cost, cost2 = dijkstra(f, m, Bc);
78   return lbl;

2. Testando as funções

 1 from result_6_twsm import tester
 2 import activity_oalexandrefelipe_6_twsm as my
 3 import ia636, ia870;
 4 import numpy as np;
 5 Bc = np.array([[0,1,0],[1,0,1], [0,1,0]]);
 6 f = adreadgray('lenina.pgm');
 7 g = ia870.iagradm(f);
 8 m = np.zeros_like(g);
 9 m[100:110, 10:20] = 1;
10 m[-20:-10, -20:-10] = 2;
11 m[20:30, -20:-10] = 3;
12 
13 
14 lbl, cost, cost2 = my.dijkstra(g, m, Bc);
15 adshow(ia636.ianormalize(g), 'Gradient image');
16 adshow(ia870.iaglblshow(m), 'Markers');
17 
18 adshow(ia636.ianormalize(cost), 'Minimum cost to one marker');
19 adshow(ia636.ianormalize(cost2 ), 'Second cost component.');
20 adshow(ia870.iaglblshow(lbl), 'Segmentation');

Gradient image

Markers

Minimum cost to one marker

Second cost component.

Segmentation

Testes numéricos, 1 casos:

 1 for i in range(1):
 2 
 3     f,m,Bc = tester.get_input(i)
 4     k1 = tester.get_output(i)
 5 
 6     print
 7     print 'Caso ', i
 8     print
 9     print 'f=\n', f * 1
10     print 'm=\n', m
11     print 'Bc=\n', Bc * 1
12     print 'Resultado esperado =\n', k1
13     g1 = my.twsm(f,m,Bc)
14     print 'Meu resultado =\n', g1
Caso  0

f=
[[10 10 10 10 10 10 10]
 [10  9  6 18  6  5 10]
 [10  9  6 18  6  8 10]
 [10  9  9 15  9  9 10]
 [10  9  9 15 12 10 10]
 [10 10 10  8 10 10 10]]
m=
[[0 0 0 0 0 0 0]
 [0 0 1 0 0 2 0]
 [0 0 1 0 0 0 0]
 [0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0]
 [0 0 0 3 0 0 0]]
Bc=
[[0 1 0]
 [1 0 1]
 [0 1 0]]
Resultado esperado =
[[1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 0 2 2]
 [1 1 0 3 3 0 0]]
Meu resultado =
[[1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 2 2 2]
 [1 1 1 0 0 2 2]
 [1 1 0 3 3 0 0]]

Testes com imagens, 1 caso:

 1 import numpy as np
 2 import ia636,ia870
 3 for i in range(1,2):
 4     print
 5     print 'Caso', i
 6     print
 7     f,m,Bc = tester.get_input(i)
 8     adshow(f, 'f.shape=%s caso %d' % (f.shape, i) )
 9     adshow(ia870.iaglblshow(m), 'marcador rotulado' )
10     print 'Bc=\n',Bc * 1
11     k1 = tester.get_output(i)
12     adshow(ia870.iaglblshow(k1), 'esperado, k1.shape=%s caso %d' % (k1.shape,i,) )
13     kmy1, c1, c2 = my.dijkstra(f,m,Bc)
14     adshow(ia870.iaglblshow(kmy1), 'meu, kmy1.shape=%s caso %d' % (kmy1.shape,i,) )
15     adshow(ia636.ianormalize(c1), 'Custo primário caso %d.' % i);
16     adshow(ia636.ianormalize(c2), 'Custo secundário caso %d.' % i);
17     adshow(ia636.ianormalize(c1[-50:,:50]), 'Custo primário caso %d.' % i);
18     adshow(ia636.ianormalize(np.where(c2[-50:,:50]> 10, 10,
19        c2[-50:,:50])), 'Custo secundário caso %d.' % i);
Caso 1

Bc=
[[1 1 1]
 [1 0 1]
 [1 1 1]]

f.shape=(128, 128) caso 1

marcador rotulado

esperado, k1.shape=(128, 128) caso 1

meu, kmy1.shape=(128, 128) caso 1

Custo primário caso 1.

Custo secundário caso 1.

Custo primário caso 1.

Custo secundário caso 1.

3. Colocando a função no sistema de testes do Adessowiki

Uma vez feita a função que você considera que esta funcionando, este trecho abaixo irá registrá-lo no teste automático, aparecendo na página resultados.

1 # Na linha abaixo, adiciona-se a função a ser testado. Caso queira registrar mais de uma
2 # função, é só adicionar mais uma linha com o run_test:
3 tester.run_test(my.twsm,'oalexandrefelipe')
4 tester.show_results_table()
Tabela de resultados consolidados de 0 funções distribuídas em 0 páginas distintas.
Ranking do autor Ranking da função Autor Função num 1 (tempo em ms) num 1 (pontuação) img 1 (tempo em ms) img 1 (pontuação) Tempo total (ms) Pontuação total
1 1 oalexandrefelipe twsm() 5.631 100% 1986.952 100% 1992.583 100%
[–] Comments
Alexandre Felipe at 2015-05-27 10:58h
 0
Obrigado professor, a sua explicação me fez entender...

Roberto Lotufo at 2015-05-26 23:29h
 0
Olá Alexandre, estes algoritmos são bastante complexos. Deve existir algum erro ainda no tratamento do segundo componente do custo. Não é fácil analisar. O mais simples é criar um exemplo numérico simples para mostrar a discrepância entre as soluções. O projeto de exemplos de teste é uma parte importante da solução.

Alexandre Felipe at 2015-05-26 23:14h
 0
Professor, preciso de ajuda de novo!

Roberto Lotufo at 2015-05-26 14:58h
 0
Olá Alexandre, a diferença que está ocorrendo é que o custo implementado no tie-zone relatado no artigo. O custo é composto de duas partes, uma é o valor do pixel (gradiente) e a outra componente do custo é quando o valor do pixel é o mesmo (platô) que é implementado intrinsicamente pela FIFO da fila hierárquica. Isto é o que torna este algoritmo do tie-zone mais complexo. Se não fosse isso, seu algoritmo estaria correto.

Alexandre Felipe at 2015-05-26 14:11h
 0
Olá professor, você tem alguma indicação sobre o que está errado agora? Eu escrevi o algoritmo de caminho mínimo mais preocupado com ter um código limpo que com desempenho. Usei a variável Bci para andar pelos vizinhos de um dado pixel, testo se o pixel está dentro da imagem, em seguida fiz o que eu entendo ser um algoritmo para busca de caminhos mínimos. Parece que meu erro está na determinação da zona de empate. Outra dúvida que tenho, como fica o final do caminho, o píxel de destino conta como parte do caminho (foi o que eu implementei), e o pixel de origem (nessa implementação ele está incluído)

Alexandre Felipe at 2015-05-26 14:06h
 0
Ok, isso deixou o resultado mais bonito para se visualizar.

Roberto Lotufo at 2015-05-26 13:04h
 0
Olá Alexandre, o custo é dado pelo maior peso no caminho, normalmente o peso é dado por uma dissimilaridade entre pixels. A forma mais fácil de calcular esta dissimilaridade é o gradiente. Assim, o watershed roda sobre a imagem do gradiente e não sobre os pixels da imagem original.

Alexandre Felipe at 2015-05-26 10:26h
 0
Sim, professor eu já tinha notado isso, eu só estava tentando bolar uma heurística para paralelizar a busca, para ver se em termos de velocidade seria melhor. Quanto a ideia do gradiente ainda não olhei para ela, estou tentando seguir diretamente pela definição. A definição que eu havia entendido é que cada ponto é declarado com a mesma classe do marcador que tenha o caminho de menor custo até o ponto. Sendo que o custo do caminho é o máximo dos valores dos pixels usados (visitados).

Roberto Lotufo at 2015-05-25 17:28h
 0
Olá Alexandre, o custo do watershed é definido como um custo simultâneo no grafo. Há necessidade de se criar uma floresta. Neste sentido sua solução de calcular os custos de forma não simultânea não funciona. Adicionalmente, como sugestão, a image do gradiente, assim, troque a imagem ilustrativa lenina.pgm pelo gradiente: gradm(f)