Filtering the pixel oriented max-tree - Wilkinson et al. Algorithm

MainPage

 1 #include<iostream>
 2 #include<vector>
 3 #include <omp.h>
 4 
 5 void direct_filter_wilkinson(int nthreads,double lamb, int h_par, int *par, int h_attr, double *attr,
 6      int h_f, unsigned char *f, int h_filt, int *filt){
 7 
 8    int TID,w, u,v;
 9    unsigned char val;
10    if (nthreads > omp_get_max_threads())
11       nthreads = omp_get_max_threads();
12 
13    omp_set_num_threads(nthreads);
14 
15    int chunk = (h_par-1)/nthreads + 1;
16    //#pragma omp parallel for shared(lamb, par, attr, f, filt) private(u,w,val)
17    #pragma omp parallel shared(lamb, h_par,par, attr, f, filt,chunk) private(u,w,val,v,TID)
18    {
19 
20       TID = omp_get_thread_num();
21       int v_min = TID*chunk;
22       int v_max = v_min + chunk;
23       if (v_max >= h_par)
24          v_max = h_par;
25 
26       for (v = v_min; v < v_max; v++){
27          if (filt[v] == -1){
28             w = v;
29             while (par[w] !=w && filt[w] == -1 && (f[w] == f[par[w]] || attr[w] < lamb)){
30                w = par[w];
31                }
32             if (filt[w]!= -1) val = filt[w];
33             else if (attr[w] >= lamb) val = f[w];
34             else val = 0;
35 
36             u = v;
37             while (u!= w){
38                if (u>= v_min && u<v_max)
39                   filt[u] = val;
40                u = par[u];
41 
42                }
43             if (w>= v_min && w<v_max)
44                filt[w] = val;
45          }
46       }
47 }
48 
49    }
OK

Interface

1 void direct_filter_wilkinson(int nthreads,double lamb, int DIM1, int *INPLACE_ARRAY1, int DIM1, double *INPLACE_ARRAY1,
2                                int DIM1, unsigned char *INPLACE_ARRAY1, int DIM1, int *INPLACE_ARRAY1);
OK [C/C++ extension is up-to-date]

Python wrapper

 1 import numpy as np
 2 
 3 from pixel_oriented_filering2 import direct_filter_wilkinson
 4 
 5 def direct_filter2(lambda_,parent,img,attr, nthreads = 1):
 6     flat_img = np.ascontiguousarray(img.ravel())
 7     out = np.empty(parent.size, dtype = np.int32)
 8     out[:] = -1
 9     direct_filter_wilkinson(nthreads,lambda_,parent,attr.astype(float),flat_img,out)
10     return out.reshape(img.shape)

Example

 1 from max_tree_alpha import MaxTreeAlpha
 2 from pixel_oriented_filering import  compute_area
 3 from pixel_oriented_filering2 import  direct_filter2
 4 import numpy as np
 5 import time
 6 from max_tree_c_01 import build_max_tree
 7 
 8 img = adreadgray('nplate/IMG_0903.jpg')
 9 #img = adreadgray('MVBook/bamboo.png')
10 #img = adreadgray('cameraman.tif')
11 Bc = np.ones((3,3), dtype = bool)
12 
13 t = time.time()
14 par,S_rev = build_max_tree(img,Bc, option = 0)
15 t = time.time() - t
16 print "MaxTree (parent,S)"
17 print "Construction time:", t
18 
19 
20 area = compute_area(S_rev,par)
21 t = time.time()
22 out = direct_filter2(200.0,par,img,area, nthreads = 2)
23 t = time.time() - t
24 print "Filtering time pixel oriented max-tree (2 threads):", t
25 
26 
27 t = time.time()
28 par,S_rev = build_max_tree(img,Bc, option = 0)
29 t = time.time() - t
30 print "MaxTree (parent,S)"
31 print "Construction time:", t
32 
33 
34 area = compute_area(S_rev,par)
35 t = time.time()
36 out = direct_filter2(200.0,par,img,area, nthreads = 4)
37 t = time.time() - t
38 print "Filtering time pixel oriented max-tree (4 threads):", t
39 
40 t = time.time()
41 par,S_rev = build_max_tree(img,Bc, option = 0)
42 t = time.time() - t
43 print "MaxTree (parent,S)"
44 print "Construction time:", t
45 
46 
47 area = compute_area(S_rev,par)
48 t = time.time()
49 out = direct_filter2(200.0,par,img,area, nthreads = 8)
50 t = time.time() - t
51 print "Filtering time pixel oriented max-tree (8 threads):", t
52 
53 t = time.time()
54 par,S_rev = build_max_tree(img,Bc, option = 0)
55 t = time.time() - t
56 print "MaxTree (parent,S)"
57 print "Construction time:", t
58 
59 
60 area = compute_area(S_rev,par)
61 t = time.time()
62 out = direct_filter2(200.0,par,img,area)
63 t = time.time() - t
64 print "Filtering time pixel oriented max-tree (1 thread):", t
65 
66 
67 
68 
69 t = time.time()
70 mxt = MaxTreeAlpha(img,Bc)
71 t = time.time() - t
72 print "MaxTreeAlpha"
73 print "Construction time:", t
74 
75 t = time.time()
76 mxt.areaOpen(200)
77 out2 = mxt.getImage()
78 t = time.time() - t
79 print "Filtering and reconstruction time:", t
80 
81 print "Os resultados são iguais:", (out == out2).all()
82 adshow(img, "Imagem original")
83 adshow(out,"Imagem filtrada (area-open)")
MaxTree (parent,S)
Construction time: 2.80575609207
Filtering time pixel oriented max-tree (2 threads): 0.268739938736
MaxTree (parent,S)
Construction time: 2.50226402283
Filtering time pixel oriented max-tree (4 threads): 0.284893989563
MaxTree (parent,S)
Construction time: 2.47825479507
Filtering time pixel oriented max-tree (8 threads): 0.338929891586
MaxTree (parent,S)
Construction time: 2.67105817795
Filtering time pixel oriented max-tree (1 thread): 0.333519935608
MaxTreeAlpha
Construction time: 4.24896001816
Filtering and reconstruction time: 0.18064904213
Os resultados são iguais: True

Imagem original

Imagem filtrada (area-open)