1 // Copyright NVIDIA Corporation 2007 -- Ignacio Castano <icastano@nvidia.com>
2 //
3 // Permission is hereby granted, free of charge, to any person
4 // obtaining a copy of this software and associated documentation
5 // files (the "Software"), to deal in the Software without
6 // restriction, including without limitation the rights to use,
7 // copy, modify, merge, publish, distribute, sublicense, and/or sell
8 // copies of the Software, and to permit persons to whom the
9 // Software is furnished to do so, subject to the following
10 // conditions:
11 //
12 // The above copyright notice and this permission notice shall be
13 // included in all copies or substantial portions of the Software.
14 //
15 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
17 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
19 // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
20 // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
21 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
22 // OTHER DEALINGS IN THE SOFTWARE.
23
24 // Math functions and operators to be used with vector types.
25
26 #ifndef CUDAMATH_H
27 #define CUDAMATH_H
28
29 #include <float.h>
30
31
32 inline __device__ __host__ float3 operator *(float3 a, float3 b)
33 {
34 return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
35 }
36
37 inline __device__ __host__ float3 operator *(float f, float3 v)
38 {
39 return make_float3(v.x*f, v.y*f, v.z*f);
40 }
41
42 inline __device__ __host__ float3 operator *(float3 v, float f)
43 {
44 return make_float3(v.x*f, v.y*f, v.z*f);
45 }
46
47 inline __device__ __host__ float3 operator +(float3 a, float3 b)
48 {
49 return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
50 }
51
52 inline __device__ __host__ void operator +=(float3 & b, float3 a)
53 {
54 b.x += a.x;
55 b.y += a.y;
56 b.z += a.z;
57 }
58
59 inline __device__ __host__ float3 operator -(float3 a, float3 b)
60 {
61 return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
62 }
63
64 inline __device__ __host__ void operator -=(float3 & b, float3 a)
65 {
66 b.x -= a.x;
67 b.y -= a.y;
68 b.z -= a.z;
69 }
70
71 inline __device__ __host__ float3 operator /(float3 v, float f)
72 {
73 float inv = 1.0f / f;
74 return v * inv;
75 }
76
77 inline __device__ __host__ void operator /=(float3 & b, float f)
78 {
79 float inv = 1.0f / f;
80 b.x *= inv;
81 b.y *= inv;
82 b.z *= inv;
83 }
84
85 inline __device__ __host__ bool operator ==(float3 a, float3 b)
86 {
87 return a.x == b.x && a.y == b.y && a.z == b.z;
88 }
89
dot(float3 a,float3 b)90 inline __device__ __host__ float dot(float3 a, float3 b)
91 {
92 return a.x * b.x + a.y * b.y + a.z * b.z;
93 }
94
dot(float4 a,float4 b)95 inline __device__ __host__ float dot(float4 a, float4 b)
96 {
97 return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
98 }
99
clamp(float f,float a,float b)100 inline __device__ __host__ float clamp(float f, float a, float b)
101 {
102 return max(a, min(f, b));
103 }
104
clamp(float3 v,float a,float b)105 inline __device__ __host__ float3 clamp(float3 v, float a, float b)
106 {
107 return make_float3(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b));
108 }
109
clamp(float3 v,float3 a,float3 b)110 inline __device__ __host__ float3 clamp(float3 v, float3 a, float3 b)
111 {
112 return make_float3(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z));
113 }
114
115
normalize(float3 v)116 inline __device__ __host__ float3 normalize(float3 v)
117 {
118 float len = 1.0f / sqrtf(dot(v, v));
119 return make_float3(v.x * len, v.y * len, v.z * len);
120 }
121
122
123
124
125 // Use power method to find the first eigenvector.
126 // http://www.miislita.com/information-retrieval-tutorial/matrix-tutorial-3-eigenvalues-eigenvectors.html
firstEigenVector(float matrix[6])127 inline __device__ __host__ float3 firstEigenVector( float matrix[6] )
128 {
129 // 8 iterations seems to be more than enough.
130
131 float3 row0 = make_float3(matrix[0], matrix[1], matrix[2]);
132 float3 row1 = make_float3(matrix[1], matrix[3], matrix[4]);
133 float3 row2 = make_float3(matrix[2], matrix[4], matrix[5]);
134
135 float r0 = dot(row0, row0);
136 float r1 = dot(row1, row1);
137 float r2 = dot(row2, row2);
138
139 float3 v;
140 if (r0 > r1 && r0 > r2) v = row0;
141 else if (r1 > r2) v = row1;
142 else v = row2;
143
144 //float3 v = make_float3(1.0f, 1.0f, 1.0f);
145 for(int i = 0; i < 8; i++) {
146 float x = v.x * matrix[0] + v.y * matrix[1] + v.z * matrix[2];
147 float y = v.x * matrix[1] + v.y * matrix[3] + v.z * matrix[4];
148 float z = v.x * matrix[2] + v.y * matrix[4] + v.z * matrix[5];
149 float m = max(max(x, y), z);
150 float iv = 1.0f / m;
151 if (m == 0.0f) iv = 0.0f;
152 v = make_float3(x*iv, y*iv, z*iv);
153 }
154
155 return v;
156 }
157
singleColor(const float3 * colors)158 inline __device__ bool singleColor(const float3 * colors)
159 {
160 #if __DEVICE_EMULATION__
161 bool sameColor = false;
162 for (int i = 0; i < 16; i++)
163 {
164 sameColor &= (colors[i] == colors[0]);
165 }
166 return sameColor;
167 #else
168 __shared__ int sameColor[16];
169
170 const int idx = threadIdx.x;
171
172 sameColor[idx] = (colors[idx] == colors[0]);
173 sameColor[idx] &= sameColor[idx^8];
174 sameColor[idx] &= sameColor[idx^4];
175 sameColor[idx] &= sameColor[idx^2];
176 sameColor[idx] &= sameColor[idx^1];
177
178 return sameColor[0];
179 #endif
180 }
181
colorSums(const float3 * colors,float3 * sums)182 inline __device__ void colorSums(const float3 * colors, float3 * sums)
183 {
184 #if __DEVICE_EMULATION__
185 float3 color_sum = make_float3(0.0f, 0.0f, 0.0f);
186 for (int i = 0; i < 16; i++)
187 {
188 color_sum += colors[i];
189 }
190
191 for (int i = 0; i < 16; i++)
192 {
193 sums[i] = color_sum;
194 }
195 #else
196
197 const int idx = threadIdx.x;
198
199 sums[idx] = colors[idx];
200 sums[idx] += sums[idx^8];
201 sums[idx] += sums[idx^4];
202 sums[idx] += sums[idx^2];
203 sums[idx] += sums[idx^1];
204
205 #endif
206 }
207
bestFitLine(const float3 * colors,float3 color_sum,float3 colorMetric)208 inline __device__ float3 bestFitLine(const float3 * colors, float3 color_sum, float3 colorMetric)
209 {
210 // Compute covariance matrix of the given colors.
211 #if __DEVICE_EMULATION__
212 float covariance[6] = {0, 0, 0, 0, 0, 0};
213 for (int i = 0; i < 16; i++)
214 {
215 float3 a = (colors[i] - color_sum * (1.0f / 16.0f)) * colorMetric;
216 covariance[0] += a.x * a.x;
217 covariance[1] += a.x * a.y;
218 covariance[2] += a.x * a.z;
219 covariance[3] += a.y * a.y;
220 covariance[4] += a.y * a.z;
221 covariance[5] += a.z * a.z;
222 }
223 #else
224
225 const int idx = threadIdx.x;
226
227 float3 diff = (colors[idx] - color_sum * (1.0f / 16.0f)) * colorMetric;
228
229 // @@ Eliminate two-way bank conflicts here.
230 // @@ It seems that doing that and unrolling the reduction doesn't help...
231 __shared__ float covariance[16*6];
232
233 covariance[6 * idx + 0] = diff.x * diff.x; // 0, 6, 12, 2, 8, 14, 4, 10, 0
234 covariance[6 * idx + 1] = diff.x * diff.y;
235 covariance[6 * idx + 2] = diff.x * diff.z;
236 covariance[6 * idx + 3] = diff.y * diff.y;
237 covariance[6 * idx + 4] = diff.y * diff.z;
238 covariance[6 * idx + 5] = diff.z * diff.z;
239
240 for(int d = 8; d > 0; d >>= 1)
241 {
242 if (idx < d)
243 {
244 covariance[6 * idx + 0] += covariance[6 * (idx+d) + 0];
245 covariance[6 * idx + 1] += covariance[6 * (idx+d) + 1];
246 covariance[6 * idx + 2] += covariance[6 * (idx+d) + 2];
247 covariance[6 * idx + 3] += covariance[6 * (idx+d) + 3];
248 covariance[6 * idx + 4] += covariance[6 * (idx+d) + 4];
249 covariance[6 * idx + 5] += covariance[6 * (idx+d) + 5];
250 }
251 }
252
253 #endif
254
255 // Compute first eigen vector.
256 return firstEigenVector(covariance);
257 }
258
259
260 #endif // CUDAMATH_H
261