Fix T46403: motion tracking not workig with Xcode 7 on OS X.
[blender.git] / extern / libmv / libmv / image / array_nd.h
1 // Copyright (c) 2007, 2008 libmv authors.
2 //
3 // Permission is hereby granted, free of charge, to any person obtaining a copy
4 // of this software and associated documentation files (the "Software"), to
5 // deal in the Software without restriction, including without limitation the
6 // rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
7 // sell copies of the Software, and to permit persons to whom the Software is
8 // furnished to do so, subject to the following conditions:
9 //
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
12 //
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
18 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
19 // IN THE SOFTWARE.
20
21 #ifndef LIBMV_IMAGE_ARRAY_ND_H
22 #define LIBMV_IMAGE_ARRAY_ND_H
23
24 #include <cassert>
25 #include <cstdio>
26 #include <cstring>
27
28 #include "libmv/image/tuple.h"
29
30 namespace libmv {
31
32 class BaseArray {};
33
34 /// A multidimensional array class.
35 template <typename T, int N>
36 class ArrayND : public BaseArray {
37  public:
38   typedef T Scalar;
39
40   /// Type for the multidimensional indices.
41   typedef Tuple<int, N> Index;
42
43   /// Create an empty array.
44   ArrayND() : data_(NULL), own_data_(true) { Resize(Index(0)); }
45
46   /// Create an array with the specified shape.
47   ArrayND(const Index &shape) : data_(NULL), own_data_(true) { Resize(shape); }
48
49   /// Create an array with the specified shape.
50   ArrayND(int *shape) : data_(NULL), own_data_(true) { Resize(shape); }
51
52   /// Copy constructor.
53   ArrayND(const ArrayND<T, N> &b) : data_(NULL), own_data_(true) {
54     ResizeLike(b);
55     std::memcpy(Data(), b.Data(), sizeof(T) * Size());
56   }
57
58   ArrayND(int s0) : data_(NULL), own_data_(true) { Resize(s0); }
59   ArrayND(int s0, int s1) : data_(NULL), own_data_(true) { Resize(s0, s1); }
60   ArrayND(int s0, int s1, int s2) : data_(NULL), own_data_(true) {
61      Resize(s0, s1, s2);
62   }
63
64   ArrayND(T* data, int s0, int s1, int s2)
65       : shape_(0), strides_(0), data_(data), own_data_(false) {
66     Resize(s0, s1, s2);
67   }
68
69   /// Destructor deletes pixel data.
70   ~ArrayND() {
71     if (own_data_) {
72       delete [] data_;
73     }
74   }
75
76   /// Assignation copies pixel data.
77   ArrayND &operator=(const ArrayND<T, N> &b) {
78     assert(this != &b);
79     ResizeLike(b);
80     std::memcpy(Data(), b.Data(), sizeof(T) * Size());
81     return *this;
82   }
83
84   const Index &Shapes() const {
85     return shape_;
86   }
87
88   const Index &Strides() const {
89     return strides_;
90   }
91
92   /// Create an array of shape s.
93   void Resize(const Index &new_shape) {
94     if (data_ != NULL && shape_ == new_shape) {
95       // Don't bother realloacting if the shapes match.
96       return;
97     }
98     shape_.Reset(new_shape);
99     strides_(N - 1) = 1;
100     for (int i = N - 1; i > 0; --i) {
101       strides_(i - 1) = strides_(i) * shape_(i);
102     }
103     if (own_data_) {
104       delete [] data_;
105       data_ = NULL;
106       if (Size() > 0) {
107         data_ = new T[Size()];
108       }
109     }
110   }
111
112   template<typename D>
113   void ResizeLike(const ArrayND<D, N> &other) {
114     Resize(other.Shape());
115   }
116
117   /// Resizes the array to shape s.  All data is lost.
118   void Resize(const int *new_shape_array) {
119     Resize(Index(new_shape_array));
120   }
121
122   /// Resize a 1D array to length s0.
123   void Resize(int s0) {
124     assert(N == 1);
125     int shape[] = {s0};
126     Resize(shape);
127   }
128
129   /// Resize a 2D array to shape (s0,s1).
130   void Resize(int s0, int s1) {
131     int shape[N] = {s0, s1};
132     for (int i = 2; i < N; ++i) {
133       shape[i] = 1;
134     }
135     Resize(shape);
136   }
137
138   // Match Eigen2's API.
139   void resize(int rows, int cols) {
140     Resize(rows, cols);
141   }
142
143   /// Resize a 3D array to shape (s0,s1,s2).
144   void Resize(int s0, int s1, int s2) {
145     assert(N == 3);
146     int shape[] = {s0, s1, s2};
147     Resize(shape);
148   }
149
150   template<typename D>
151   void CopyFrom(const ArrayND<D, N> &other) {
152     ResizeLike(other);
153     T *data = Data();
154     const D *other_data = other.Data();
155     for (int i = 0; i < Size(); ++i) {
156       data[i] = T(other_data[i]);
157     }
158   }
159
160   void Fill(T value) {
161     for (int i = 0; i < Size(); ++i) {
162       Data()[i] = value;
163     }
164   }
165
166   // Match Eigen's API.
167   void fill(T value) {
168     for (int i = 0; i < Size(); ++i) {
169       Data()[i] = value;
170     }
171   }
172
173   /// Return a tuple containing the length of each axis.
174   const Index &Shape() const {
175     return shape_;
176   }
177
178   /// Return the length of an axis.
179   int Shape(int axis) const {
180     return shape_(axis);
181   }
182
183   /// Return the distance between neighboring elements along axis.
184   int Stride(int axis) const {
185     return strides_(axis);
186   }
187
188   /// Return the number of elements of the array.
189   int Size() const {
190     int size = 1;
191     for (int i = 0; i < N; ++i)
192       size *= Shape(i);
193     return size;
194   }
195
196   /// Return the total amount of memory used by the array.
197   int MemorySizeInBytes() const {
198     return sizeof(*this) + Size() * sizeof(T);
199   }
200
201   /// Pointer to the first element of the array.
202   T *Data() { return data_; }
203
204   /// Constant pointer to the first element of the array.
205   const T *Data() const { return data_; }
206
207   /// Distance between the first element and the element at position index.
208   int Offset(const Index &index) const {
209     int offset = 0;
210     for (int i = 0; i < N; ++i)
211       offset += index(i) * Stride(i);
212     return offset;
213   }
214
215   /// 1D specialization.
216   int Offset(int i0) const {
217     assert(N == 1);
218     return i0 * Stride(0);
219   }
220
221   /// 2D specialization.
222   int Offset(int i0, int i1) const {
223     assert(N == 2);
224     return i0 * Stride(0) + i1 * Stride(1);
225   }
226
227   /// 3D specialization.
228   int Offset(int i0, int i1, int i2) const {
229     assert(N == 3);
230     return i0 * Stride(0) + i1 * Stride(1) + i2 * Stride(2);
231   }
232
233   /// Return a reference to the element at position index.
234   T &operator()(const Index &index) {
235     // TODO(pau) Boundary checking in debug mode.
236     return *( Data() + Offset(index) );
237   }
238
239   /// 1D specialization.
240   T &operator()(int i0) {
241     return *( Data() + Offset(i0) );
242   }
243
244   /// 2D specialization.
245   T &operator()(int i0, int i1) {
246     assert(0 <= i0 && i0 < Shape(0));
247     assert(0 <= i1 && i1 < Shape(1));
248     return *(Data() + Offset(i0, i1));
249   }
250
251   /// 3D specialization.
252   T &operator()(int i0, int i1, int i2) {
253     assert(0 <= i0 && i0 < Shape(0));
254     assert(0 <= i1 && i1 < Shape(1));
255     assert(0 <= i2 && i2 < Shape(2));
256     return *(Data() + Offset(i0, i1, i2));
257   }
258
259   /// Return a constant reference to the element at position index.
260   const T &operator()(const Index &index) const {
261     return *(Data() + Offset(index));
262   }
263
264   /// 1D specialization.
265   const T &operator()(int i0) const {
266     return *(Data() + Offset(i0));
267   }
268
269   /// 2D specialization.
270   const T &operator()(int i0, int i1) const {
271     assert(0 <= i0 && i0 < Shape(0));
272     assert(0 <= i1 && i1 < Shape(1));
273     return *(Data() + Offset(i0, i1));
274   }
275
276   /// 3D specialization.
277   const T &operator()(int i0, int i1, int i2) const {
278     return *(Data() + Offset(i0, i1, i2));
279   }
280
281   /// True if index is inside array.
282   bool Contains(const Index &index) const {
283     for (int i = 0; i < N; ++i)
284       if (index(i) < 0 || index(i) >= Shape(i))
285         return false;
286     return true;
287   }
288
289   /// 1D specialization.
290   bool Contains(int i0) const {
291     return 0 <= i0 && i0 < Shape(0);
292   }
293
294   /// 2D specialization.
295   bool Contains(int i0, int i1) const {
296     return 0 <= i0 && i0 < Shape(0)
297         && 0 <= i1 && i1 < Shape(1);
298   }
299
300   /// 3D specialization.
301   bool Contains(int i0, int i1, int i2) const {
302     return 0 <= i0 && i0 < Shape(0)
303         && 0 <= i1 && i1 < Shape(1)
304         && 0 <= i2 && i2 < Shape(2);
305   }
306
307   bool operator==(const ArrayND<T, N> &other) const {
308     if (shape_ != other.shape_) return false;
309     if (strides_ != other.strides_) return false;
310     for (int i = 0; i < Size(); ++i) {
311       if (this->Data()[i] != other.Data()[i])
312         return false;
313     }
314     return true;
315   }
316
317   bool operator!=(const ArrayND<T, N> &other) const {
318     return !(*this == other);
319   }
320
321   ArrayND<T, N> operator*(const ArrayND<T, N> &other) const {
322     assert(Shape() = other.Shape());
323     ArrayND<T, N> res;
324     res.ResizeLike(*this);
325     for (int i = 0; i < res.Size(); ++i) {
326       res.Data()[i] = Data()[i] * other.Data()[i];
327     }
328     return res;
329   }
330
331  protected:
332   /// The number of element in each dimension.
333   Index shape_;
334
335   /// How to jump to neighbors in each dimension.
336   Index strides_;
337
338   /// Pointer to the first element of the array.
339   T *data_;
340
341   /// Flag if this Array either own or reference the data
342   bool own_data_;
343 };
344
345 /// 3D array (row, column, channel).
346 template <typename T>
347 class Array3D : public ArrayND<T, 3> {
348   typedef ArrayND<T, 3> Base;
349  public:
350   Array3D()
351       : Base() {
352   }
353   Array3D(int height, int width, int depth = 1)
354       : Base(height, width, depth) {
355   }
356   Array3D(T* data, int height, int width, int depth = 1)
357       : Base(data, height, width, depth) {
358   }
359
360   void Resize(int height, int width, int depth = 1) {
361     Base::Resize(height, width, depth);
362   }
363
364   int Height() const {
365     return Base::Shape(0);
366   }
367   int Width() const {
368     return Base::Shape(1);
369   }
370   int Depth() const {
371     return Base::Shape(2);
372   }
373
374   // Match Eigen2's API so that Array3D's and Mat*'s can work together via
375   // template magic.
376   int rows() const { return Height(); }
377   int cols() const { return Width(); }
378   int depth() const { return Depth(); }
379
380   int Get_Step() const { return Width()*Depth(); }
381
382   /// Enable accessing with 2 indices for grayscale images.
383   T &operator()(int i0, int i1, int i2 = 0) {
384     assert(0 <= i0 && i0 < Height());
385     assert(0 <= i1 && i1 < Width());
386     return Base::operator()(i0, i1, i2);
387   }
388   const T &operator()(int i0, int i1, int i2 = 0) const {
389     assert(0 <= i0 && i0 < Height());
390     assert(0 <= i1 && i1 < Width());
391     return Base::operator()(i0, i1, i2);
392   }
393 };
394
395 typedef Array3D<unsigned char> Array3Du;
396 typedef Array3D<unsigned int> Array3Dui;
397 typedef Array3D<int> Array3Di;
398 typedef Array3D<float> Array3Df;
399 typedef Array3D<short> Array3Ds;
400
401 void SplitChannels(const Array3Df &input,
402                    Array3Df *channel0,
403                    Array3Df *channel1,
404                    Array3Df *channel2);
405
406 void PrintArray(const Array3Df &array);
407
408 /** Convert a float array into a byte array by scaling values by 255* (max-min).
409  *  where max and min are automatically detected 
410  *  (if automatic_range_detection = true)
411  * \note and TODO this automatic detection only works when the image contains
412  *  at least one pixel of both bounds.
413  **/
414 void FloatArrayToScaledByteArray(const Array3Df &float_array,
415                                  Array3Du *byte_array,
416                                  bool automatic_range_detection = false);
417
418 //! Convert a byte array into a float array by dividing values by 255.
419 void ByteArrayToScaledFloatArray(const Array3Du &byte_array,
420                                  Array3Df *float_array);
421
422 template <typename AArrayType, typename BArrayType, typename CArrayType>
423 void MultiplyElements(const AArrayType &a,
424            const BArrayType &b,
425            CArrayType *c) {
426   // This function does an element-wise multiply between
427   // the two Arrays A and B, and stores the result in C.
428   // A and B must have the same dimensions.
429   assert(a.Shape() == b.Shape());
430   c->ResizeLike(a);
431
432   // To perform the multiplcation, a "current" index into the N-dimensions of
433   // the A and B matrix specifies which elements are being multiplied.
434   typename CArrayType::Index index;
435
436   // The index starts at the maximum value for each dimension
437   const typename CArrayType::Index& cShape = c->Shape();
438   for ( int i = 0; i < CArrayType::Index::SIZE; ++i )
439     index(i) = cShape(i) - 1;
440
441   // After each multiplication, the highest-dimensional index is reduced.
442   // if this reduces it less than zero, it resets to its maximum value
443   // and decrements the index of the next lower dimension.
444   // This ripple-action continues until the entire new array has been
445   // calculated, indicated by dimension zero having a negative index.
446   while ( index(0) >= 0 ) {
447     (*c)(index) = a(index) * b(index);
448
449     int dimension = CArrayType::Index::SIZE - 1;
450     index(dimension) = index(dimension) - 1;
451     while ( dimension > 0 && index(dimension) < 0 ) {
452       index(dimension) = cShape(dimension) - 1;
453       index(dimension - 1) = index(dimension - 1) - 1;
454       --dimension;
455     }
456   }
457 }
458
459 template <typename TA, typename TB, typename TC>
460 void MultiplyElements(const ArrayND<TA, 3> &a,
461                       const ArrayND<TB, 3> &b,
462                       ArrayND<TC, 3> *c) {
463   // Specialization for N==3
464   c->ResizeLike(a);
465   assert(a.Shape(0) == b.Shape(0));
466   assert(a.Shape(1) == b.Shape(1));
467   assert(a.Shape(2) == b.Shape(2));
468   for (int i = 0; i < a.Shape(0); ++i) {
469     for (int j = 0; j < a.Shape(1); ++j) {
470       for (int k = 0; k < a.Shape(2); ++k) {
471         (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k));
472       }
473     }
474   }
475 }
476
477 template <typename TA, typename TB, typename TC>
478 void MultiplyElements(const Array3D<TA> &a,
479                       const Array3D<TB> &b,
480                       Array3D<TC> *c) {
481   // Specialization for N==3
482   c->ResizeLike(a);
483   assert(a.Shape(0) == b.Shape(0));
484   assert(a.Shape(1) == b.Shape(1));
485   assert(a.Shape(2) == b.Shape(2));
486   for (int i = 0; i < a.Shape(0); ++i) {
487     for (int j = 0; j < a.Shape(1); ++j) {
488       for (int k = 0; k < a.Shape(2); ++k) {
489         (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k));
490       }
491     }
492   }
493 }
494
495 }  // namespace libmv
496
497 #endif  // LIBMV_IMAGE_ARRAY_ND_H