1 //---------------------------------------------------------------------------//
2 // Copyright (c) 2013-2014 Kyle Lutz <kyle.r.lutz@gmail.com>
3 //
4 // Distributed under the Boost Software License, Version 1.0
5 // See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt
7 //
8 // See http://boostorg.github.com/compute for more information.
9 //---------------------------------------------------------------------------//
10
11 #include <opencv2/core/core.hpp>
12 #include <opencv2/highgui/highgui.hpp>
13 #include <opencv2/imgproc/imgproc.hpp>
14
15 #include <boost/compute/system.hpp>
16 #include <boost/compute/container/vector.hpp>
17 #include <boost/compute/image/image2d.hpp>
18 #include <boost/compute/interop/opencv/core.hpp>
19 #include <boost/compute/interop/opencv/highgui.hpp>
20 #include <boost/compute/random/default_random_engine.hpp>
21 #include <boost/compute/random/uniform_real_distribution.hpp>
22 #include <boost/compute/utility/dim.hpp>
23 #include <boost/compute/utility/source.hpp>
24
25 namespace compute = boost::compute;
26
27 using compute::dim;
28 using compute::int_;
29 using compute::float_;
30 using compute::float2_;
31
32 // the k-means example implements the k-means clustering algorithm
main()33 int main()
34 {
35 // number of clusters
36 size_t k = 6;
37
38 // number of points
39 size_t n_points = 4500;
40
41 // height and width of image
42 size_t height = 800;
43 size_t width = 800;
44
45 // get default device and setup context
46 compute::device gpu = compute::system::default_device();
47 compute::context context(gpu);
48 compute::command_queue queue(context, gpu);
49
50 // generate random, uniformily-distributed points
51 compute::default_random_engine random_engine(queue);
52 compute::uniform_real_distribution<float_> uniform_distribution(0, 800);
53
54 compute::vector<float2_> points(n_points, context);
55 uniform_distribution.generate(
56 compute::make_buffer_iterator<float_>(points.get_buffer(), 0),
57 compute::make_buffer_iterator<float_>(points.get_buffer(), n_points * 2),
58 random_engine,
59 queue
60 );
61
62 // initialize all points to cluster 0
63 compute::vector<int_> clusters(n_points, context);
64 compute::fill(clusters.begin(), clusters.end(), 0, queue);
65
66 // create initial means with the first k points
67 compute::vector<float2_> means(k, context);
68 compute::copy_n(points.begin(), k, means.begin(), queue);
69
70 // k-means clustering program source
71 const char k_means_source[] = BOOST_COMPUTE_STRINGIZE_SOURCE(
72 __kernel void assign_clusters(__global const float2 *points,
73 __global const float2 *means,
74 const int k,
75 __global int *clusters)
76 {
77 const uint gid = get_global_id(0);
78
79 const float2 point = points[gid];
80
81 // find the closest cluster
82 float current_distance = 0;
83 int closest_cluster = -1;
84
85 // find closest cluster mean to the point
86 for(int i = 0; i < k; i++){
87 const float2 mean = means[i];
88
89 int distance_to_mean = distance(point, mean);
90 if(closest_cluster == -1 || distance_to_mean < current_distance){
91 current_distance = distance_to_mean;
92 closest_cluster = i;
93 }
94 }
95
96 // write new cluster
97 clusters[gid] = closest_cluster;
98 }
99
100 __kernel void update_means(__global const float2 *points,
101 const uint n_points,
102 __global float2 *means,
103 __global const int *clusters)
104 {
105 const uint k = get_global_id(0);
106
107 float2 sum = { 0, 0 };
108 float count = 0;
109 for(uint i = 0; i < n_points; i++){
110 if(clusters[i] == k){
111 sum += points[i];
112 count += 1;
113 }
114 }
115
116 means[k] = sum / count;
117 }
118 );
119
120 // build the k-means program
121 compute::program k_means_program =
122 compute::program::build_with_source(k_means_source, context);
123
124 // setup the k-means kernels
125 compute::kernel assign_clusters_kernel(k_means_program, "assign_clusters");
126 assign_clusters_kernel.set_arg(0, points);
127 assign_clusters_kernel.set_arg(1, means);
128 assign_clusters_kernel.set_arg(2, int_(k));
129 assign_clusters_kernel.set_arg(3, clusters);
130
131 compute::kernel update_means_kernel(k_means_program, "update_means");
132 update_means_kernel.set_arg(0, points);
133 update_means_kernel.set_arg(1, int_(n_points));
134 update_means_kernel.set_arg(2, means);
135 update_means_kernel.set_arg(3, clusters);
136
137 // run the k-means algorithm
138 for(int iteration = 0; iteration < 25; iteration++){
139 queue.enqueue_1d_range_kernel(assign_clusters_kernel, 0, n_points, 0);
140 queue.enqueue_1d_range_kernel(update_means_kernel, 0, k, 0);
141 }
142
143 // create output image
144 compute::image2d image(
145 context, width, height, compute::image_format(CL_RGBA, CL_UNSIGNED_INT8)
146 );
147
148 // program with two kernels, one to fill the image with white, and then
149 // one the draw to points calculated in coordinates on the image
150 const char draw_walk_source[] = BOOST_COMPUTE_STRINGIZE_SOURCE(
151 __kernel void draw_points(__global const float2 *points,
152 __global const int *clusters,
153 __write_only image2d_t image)
154 {
155 const uint i = get_global_id(0);
156 const float2 coord = points[i];
157
158 // map cluster number to color
159 uint4 color = { 0, 0, 0, 0 };
160 switch(clusters[i]){
161 case 0:
162 color = (uint4)(255, 0, 0, 255);
163 break;
164 case 1:
165 color = (uint4)(0, 255, 0, 255);
166 break;
167 case 2:
168 color = (uint4)(0, 0, 255, 255);
169 break;
170 case 3:
171 color = (uint4)(255, 255, 0, 255);
172 break;
173 case 4:
174 color = (uint4)(255, 0, 255, 255);
175 break;
176 case 5:
177 color = (uint4)(0, 255, 255, 255);
178 break;
179 }
180
181 // draw a 3x3 pixel point
182 for(int x = -1; x <= 1; x++){
183 for(int y = -1; y <= 1; y++){
184 if(coord.x + x > 0 && coord.x + x < get_image_width(image) &&
185 coord.y + y > 0 && coord.y + y < get_image_height(image)){
186 write_imageui(image, (int2)(coord.x, coord.y) + (int2)(x, y), color);
187 }
188 }
189 }
190 }
191
192 __kernel void fill_gray(__write_only image2d_t image)
193 {
194 const int2 coord = { get_global_id(0), get_global_id(1) };
195
196 if(coord.x < get_image_width(image) && coord.y < get_image_height(image)){
197 uint4 gray = { 15, 15, 15, 15 };
198 write_imageui(image, coord, gray);
199 }
200 }
201 );
202
203 // build the program
204 compute::program draw_program =
205 compute::program::build_with_source(draw_walk_source, context);
206
207 // fill image with dark gray
208 compute::kernel fill_kernel(draw_program, "fill_gray");
209 fill_kernel.set_arg(0, image);
210
211 queue.enqueue_nd_range_kernel(
212 fill_kernel, dim(0, 0), dim(width, height), dim(1, 1)
213 );
214
215 // draw points colored according to cluster
216 compute::kernel draw_kernel(draw_program, "draw_points");
217 draw_kernel.set_arg(0, points);
218 draw_kernel.set_arg(1, clusters);
219 draw_kernel.set_arg(2, image);
220 queue.enqueue_1d_range_kernel(draw_kernel, 0, n_points, 0);
221
222 // show image
223 compute::opencv_imshow("k-means", image, queue);
224
225 // wait and return
226 cv::waitKey(0);
227
228 return 0;
229 }
230