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
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()
```MaxTree (parent,S)