Filtering the Pixel Oriented Max-tree - Carlinet and Géraud Algorithm

MainPage

 1 #include<iostream>
 2  #include<vector>
 3  #include<new>
 4 
 5  void compute_area_c(int h_S, int *S, int h_parent, int *parent, int h_area, int *area){
 6      int p,q,p_root;
 7      p_root = S[h_S-1];
 8 
 9      for (int i = 0; i< h_S - 1; i++){
10          p = S[i];
11          if (p != p_root){
12              q = parent[p];
13              area[q]+= area[p];
14              }
15          }
16      }
17 
18 
19  void compute_height_c2(int h_S, int *S, int h_parent, int *parent,int h_img, unsigned short *flat_img,
20                        int h_height, int *height){
21     int p,aux, aux2;
22     for(int i = h_S - 1; i>= 0; i--){
23        p = S[i];
24        aux = parent[p];
25        aux2 = height[p] + flat_img[p] - flat_img[aux];
26        if (aux2 > height[aux])
27           height[aux] = aux2;
28 
29     }
30 }
31 
32 
33  void direct_filter_c(double lamb, int h_S, int *S, int h_parent, int *parent, int h_img,
34   unsigned short *flat_img, int h_out, unsigned short *out, int h_attr, double *attr ){
35      int p,q,p_root;
36      p_root = S[h_S-1];
37      if (attr[p_root] < lamb)
38          out[p_root] = 0;
39      else
40          out[p_root] = flat_img[p_root];
41      for (int i = h_S - 1; i >= 0; i--){
42          p = S[i];
43          q = parent[p];
44          if (flat_img[q] == flat_img[p])
45              out[p] = out[q];
46          else if (attr[p] < lamb)
47              out[p] = out[q];
48          else
49              out[p] = flat_img[p];
50      }
51  }
OK

Interface

1 void compute_area_c(int DIM1, int *INPLACE_ARRAY1, int DIM1, int *INPLACE_ARRAY1, int DIM1, int *INPLACE_ARRAY1);
2 void compute_height_c2(int DIM1, int *INPLACE_ARRAY1, int DIM1, int *INPLACE_ARRAY1, int DIM1, unsigned short *INPLACE_ARRAY1,
3                       int DIM1, int *INPLACE_ARRAY1);
4 void direct_filter_c(double lambda, int DIM1, int *INPLACE_ARRAY1, int DIM1, int *INPLACE_ARRAY1,
5  int DIM1, unsigned short *INPLACE_ARRAY1, int DIM1, unsigned short *INPLACE_ARRAY1, int DIM1, double *INPLACE_ARRAY1 );
OK [C/C++ extension is up-to-date]

Python wrapper

 1 import numpy as np
 2 
 3 def compute_area(S,parent):
 4     area = np.ones(S.size, dtype = np.int32)
 5     compute_area_c(S,parent,area)
 6     return area
 7 
 8 def compute_height(S,parent,f):
 9     height = np.zeros(S.size, dtype = np.int32)
10     ftype = f.dtype
11     if (ftype == np.uint8):
12         fu16 = f.astype(np.uint16)
13     else:
14         fu16 = f
15     compute_height_c2(S,parent,fu16.ravel(),height)
16     return height
17 
18 def direct_filter(lambda_,S,parent,f,attr):
19     ftype = f.dtype
20     if (ftype == np.uint8):
21         fu16 = f.astype(np.uint16)
22     else:
23         fu16 = f
24     flat_img = np.ascontiguousarray(fu16.ravel())
25     out = np.empty(S.size, dtype = np.uint16)
26     direct_filter_c(lambda_,S,parent,flat_img,out,attr.astype(float))
27     return out.reshape(f.shape).astype(ftype)

Example

 1 from max_tree_alpha import MaxTreeAlpha
 2 from pixel_oriented_filering import direct_filter, compute_area, compute_height
 3 import numpy as np
 4 import time
 5 from max_tree_c_01 import build_max_tree
 6 
 7 img = adreadgray('nplate/IMG_0903.jpg')[::2,::2]
 8 Bc = np.ones((3,3), dtype = bool)
 9 
10 t = time.time()
11 par,S_rev = build_max_tree(img,Bc, option = 0)
12 t = time.time() - t
13 print "MaxTree (parent,S)"
14 print "Construction time:", t
15 S = S_rev[::-1].copy()
16 
17 t = time.time()
18 height = compute_height(S,par,img)
19 out = direct_filter(10.1,S_rev,par,img,height)
20 t = time.time() - t
21 print "Filtering time pixel oriented max-tree:", t
22 
23 
24 
25 t = time.time()
26 mxt = MaxTreeAlpha(img,Bc)
27 t = time.time() - t
28 print "MaxTreeAlpha"
29 print "Construction time:", t
30 
31 t = time.time()
32 height = mxt.computeHeight()
33 mxt.contractDR(height > 11)
34 out2 = mxt.getImage()
35 t = time.time() - t
36 print "Filtering and reconstruction time:", t
37 
38 print "Os resultados são iguais:", (out == out2).all()
39 adshow(img, "Imagem original")
40 adshow(out,"Imagem filtrada (hmax)")
41 adshow(out2,"Imagem filtrada (hmax)")
MaxTree (parent,S)
Construction time: 0.518689155579
Filtering time pixel oriented max-tree: 0.10297203064
MaxTreeAlpha
Construction time: 0.6921210289
Filtering and reconstruction time: 0.0338439941406
Os resultados são iguais: True

Imagem original

Imagem filtrada (hmax)

Imagem filtrada (hmax)