b56a765223b54933827d55952cf7372957fe944c
[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) : data_(data), own_data_(false) {
65     Resize(s0, s1, s2);
66   }
67
68   /// Destructor deletes pixel data.
69   ~ArrayND() {
70     if (own_data_) {
71       delete [] data_;
72     }
73   }
74
75   /// Assignation copies pixel data.
76   ArrayND &operator=(const ArrayND<T, N> &b) {
77     assert(this != &b);
78     ResizeLike(b);
79     std::memcpy(Data(), b.Data(), sizeof(T) * Size());
80     return *this;
81   }
82
83   const Index &Shapes() const {
84     return shape_;
85   }
86
87   const Index &Strides() const {
88     return strides_;
89   }
90
91   /// Create an array of shape s.
92   void Resize(const Index &new_shape) {
93     if (data_ != NULL && shape_ == new_shape) {
94       // Don't bother realloacting if the shapes match.
95       return;
96     }
97     shape_.Reset(new_shape);
98     strides_(N - 1) = 1;
99     for (int i = N - 1; i > 0; --i) {
100       strides_(i - 1) = strides_(i) * shape_(i);
101     }
102     if (own_data_) {
103       delete [] data_;
104       data_ = NULL;
105       if (Size() > 0) {
106         data_ = new T[Size()];
107       }
108     }
109   }
110
111   template<typename D>
112   void ResizeLike(const ArrayND<D, N> &other) {
113     Resize(other.Shape());
114   }
115
116   /// Resizes the array to shape s.  All data is lost.
117   void Resize(const int *new_shape_array) {
118     Resize(Index(new_shape_array));
119   }
120
121   /// Resize a 1D array to length s0.
122   void Resize(int s0) {
123     assert(N == 1);
124     int shape[] = {s0};
125     Resize(shape);
126   }
127
128   /// Resize a 2D array to shape (s0,s1).
129   void Resize(int s0, int s1) {
130     int shape[N] = {s0, s1};
131     for (int i = 2; i < N; ++i) {
132       shape[i] = 1;
133     }
134     Resize(shape);
135   }
136
137   // Match Eigen2's API.
138   void resize(int rows, int cols) {
139     Resize(rows, cols);
140   }
141
142   /// Resize a 3D array to shape (s0,s1,s2).
143   void Resize(int s0, int s1, int s2) {
144     assert(N == 3);
145     int shape[] = {s0, s1, s2};
146     Resize(shape);
147   }
148
149   template<typename D>
150   void CopyFrom(const ArrayND<D, N> &other) {
151     ResizeLike(other);
152     T *data = Data();
153     const D *other_data = other.Data();
154     for (int i = 0; i < Size(); ++i) {
155       data[i] = T(other_data[i]);
156     }
157   }
158
159   void Fill(T value) {
160     for (int i = 0; i < Size(); ++i) {
161       Data()[i] = value;
162     }
163   }
164
165   // Match Eigen's API.
166   void fill(T value) {
167     for (int i = 0; i < Size(); ++i) {
168       Data()[i] = value;
169     }
170   }
171
172   /// Return a tuple containing the length of each axis.
173   const Index &Shape() const {
174     return shape_;
175   }
176
177   /// Return the length of an axis.
178   int Shape(int axis) const {
179     return shape_(axis);
180   }
181
182   /// Return the distance between neighboring elements along axis.
183   int Stride(int axis) const {
184     return strides_(axis);
185   }
186
187   /// Return the number of elements of the array.
188   int Size() const {
189     int size = 1;
190     for (int i = 0; i < N; ++i)
191       size *= Shape(i);
192     return size;
193   }
194
195   /// Return the total amount of memory used by the array.
196   int MemorySizeInBytes() const {
197     return sizeof(*this) + Size() * sizeof(T);
198   }
199
200   /// Pointer to the first element of the array.
201   T *Data() { return data_; }
202
203   /// Constant pointer to the first element of the array.
204   const T *Data() const { return data_; }
205
206   /// Distance between the first element and the element at position index.
207   int Offset(const Index &index) const {
208     int offset = 0;
209     for (int i = 0; i < N; ++i)
210       offset += index(i) * Stride(i);
211     return offset;
212   }
213
214   /// 1D specialization.
215   int Offset(int i0) const {
216     assert(N == 1);
217     return i0 * Stride(0);
218   }
219
220   /// 2D specialization.
221   int Offset(int i0, int i1) const {
222     assert(N == 2);
223     return i0 * Stride(0) + i1 * Stride(1);
224   }
225
226   /// 3D specialization.
227   int Offset(int i0, int i1, int i2) const {
228     assert(N == 3);
229     return i0 * Stride(0) + i1 * Stride(1) + i2 * Stride(2);
230   }
231
232   /// Return a reference to the element at position index.
233   T &operator()(const Index &index) {
234     // TODO(pau) Boundary checking in debug mode.
235     return *( Data() + Offset(index) );
236   }
237
238   /// 1D specialization.
239   T &operator()(int i0) {
240     return *( Data() + Offset(i0) );
241   }
242
243   /// 2D specialization.
244   T &operator()(int i0, int i1) {
245     assert(0 <= i0 && i0 < Shape(0));
246     assert(0 <= i1 && i1 < Shape(1));
247     return *(Data() + Offset(i0, i1));
248   }
249
250   /// 3D specialization.
251   T &operator()(int i0, int i1, int i2) {
252     assert(0 <= i0 && i0 < Shape(0));
253     assert(0 <= i1 && i1 < Shape(1));
254     assert(0 <= i2 && i2 < Shape(2));
255     return *(Data() + Offset(i0, i1, i2));
256   }
257
258   /// Return a constant reference to the element at position index.
259   const T &operator()(const Index &index) const {
260     return *(Data() + Offset(index));
261   }
262
263   /// 1D specialization.
264   const T &operator()(int i0) const {
265     return *(Data() + Offset(i0));
266   }
267
268   /// 2D specialization.
269   const T &operator()(int i0, int i1) const {
270     assert(0 <= i0 && i0 < Shape(0));
271     assert(0 <= i1 && i1 < Shape(1));
272     return *(Data() + Offset(i0, i1));
273   }
274
275   /// 3D specialization.
276   const T &operator()(int i0, int i1, int i2) const {
277     return *(Data() + Offset(i0, i1, i2));
278   }
279
280   /// True if index is inside array.
281   bool Contains(const Index &index) const {
282     for (int i = 0; i < N; ++i)
283       if (index(i) < 0 || index(i) >= Shape(i))
284         return false;
285     return true;
286   }
287
288   /// 1D specialization.
289   bool Contains(int i0) const {
290     return 0 <= i0 && i0 < Shape(0);
291   }
292
293   /// 2D specialization.
294   bool Contains(int i0, int i1) const {
295     return 0 <= i0 && i0 < Shape(0)
296         && 0 <= i1 && i1 < Shape(1);
297   }
298
299   /// 3D specialization.
300   bool Contains(int i0, int i1, int i2) const {
301     return 0 <= i0 && i0 < Shape(0)
302         && 0 <= i1 && i1 < Shape(1)
303         && 0 <= i2 && i2 < Shape(2);
304   }
305
306   bool operator==(const ArrayND<T, N> &other) const {
307     if (shape_ != other.shape_) return false;
308     if (strides_ != other.strides_) return false;
309     for (int i = 0; i < Size(); ++i) {
310       if (this->Data()[i] != other.Data()[i])
311         return false;
312     }
313     return true;
314   }
315
316   bool operator!=(const ArrayND<T, N> &other) const {
317     return !(*this == other);
318   }
319
320   ArrayND<T, N> operator*(const ArrayND<T, N> &other) const {
321     assert(Shape() = other.Shape());
322     ArrayND<T, N> res;
323     res.ResizeLike(*this);
324     for (int i = 0; i < res.Size(); ++i) {
325       res.Data()[i] = Data()[i] * other.Data()[i];
326     }
327     return res;
328   }
329
330  protected:
331   /// The number of element in each dimension.
332   Index shape_;
333
334   /// How to jump to neighbors in each dimension.
335   Index strides_;
336
337   /// Pointer to the first element of the array.
338   T *data_;
339
340   /// Flag if this Array either own or reference the data
341   bool own_data_;
342 };
343
344 /// 3D array (row, column, channel).
345 template <typename T>
346 class Array3D : public ArrayND<T, 3> {
347   typedef ArrayND<T, 3> Base;
348  public:
349   Array3D()
350       : Base() {
351   }
352   Array3D(int height, int width, int depth = 1)
353       : Base(height, width, depth) {
354   }
355   Array3D(T* data, int height, int width, int depth = 1)
356       : Base(data, height, width, depth) {
357   }
358
359   void Resize(int height, int width, int depth = 1) {
360     Base::Resize(height, width, depth);
361   }
362
363   int Height() const {
364     return Base::Shape(0);
365   }
366   int Width() const {
367     return Base::Shape(1);
368   }
369   int Depth() const {
370     return Base::Shape(2);
371   }
372
373   // Match Eigen2's API so that Array3D's and Mat*'s can work together via
374   // template magic.
375   int rows() const { return Height(); }
376   int cols() const { return Width(); }
377   int depth() const { return Depth(); }
378
379   int Get_Step() const { return Width()*Depth(); }
380
381   /// Enable accessing with 2 indices for grayscale images.
382   T &operator()(int i0, int i1, int i2 = 0) {
383     assert(0 <= i0 && i0 < Height());
384     assert(0 <= i1 && i1 < Width());
385     return Base::operator()(i0, i1, i2);
386   }
387   const T &operator()(int i0, int i1, int i2 = 0) const {
388     assert(0 <= i0 && i0 < Height());
389     assert(0 <= i1 && i1 < Width());
390     return Base::operator()(i0, i1, i2);
391   }
392 };
393
394 typedef Array3D<unsigned char> Array3Du;
395 typedef Array3D<unsigned int> Array3Dui;
396 typedef Array3D<int> Array3Di;
397 typedef Array3D<float> Array3Df;
398 typedef Array3D<short> Array3Ds;
399
400 void SplitChannels(const Array3Df &input,
401                    Array3Df *channel0,
402                    Array3Df *channel1,
403                    Array3Df *channel2);
404
405 void PrintArray(const Array3Df &array);
406
407 /** Convert a float array into a byte array by scaling values by 255* (max-min).
408  *  where max and min are automatically detected 
409  *  (if automatic_range_detection = true)
410  * \note and TODO this automatic detection only works when the image contains
411  *  at least one pixel of both bounds.
412  **/
413 void FloatArrayToScaledByteArray(const Array3Df &float_array,
414                                  Array3Du *byte_array,
415                                  bool automatic_range_detection = false);
416
417 //! Convert a byte array into a float array by dividing values by 255.
418 void ByteArrayToScaledFloatArray(const Array3Du &byte_array,
419                                  Array3Df *float_array);
420
421 template <typename AArrayType, typename BArrayType, typename CArrayType>
422 void MultiplyElements(const AArrayType &a,
423            const BArrayType &b,
424            CArrayType *c) {
425   // This function does an element-wise multiply between
426   // the two Arrays A and B, and stores the result in C.
427   // A and B must have the same dimensions.
428   assert(a.Shape() == b.Shape());
429   c->ResizeLike(a);
430
431   // To perform the multiplcation, a "current" index into the N-dimensions of
432   // the A and B matrix specifies which elements are being multiplied.
433   typename CArrayType::Index index;
434
435   // The index starts at the maximum value for each dimension
436   const typename CArrayType::Index& cShape = c->Shape();
437   for ( int i = 0; i < CArrayType::Index::SIZE; ++i )
438     index(i) = cShape(i) - 1;
439
440   // After each multiplication, the highest-dimensional index is reduced.
441   // if this reduces it less than zero, it resets to its maximum value
442   // and decrements the index of the next lower dimension.
443   // This ripple-action continues until the entire new array has been
444   // calculated, indicated by dimension zero having a negative index.
445   while ( index(0) >= 0 ) {
446     (*c)(index) = a(index) * b(index);
447
448     int dimension = CArrayType::Index::SIZE - 1;
449     index(dimension) = index(dimension) - 1;
450     while ( dimension > 0 && index(dimension) < 0 ) {
451       index(dimension) = cShape(dimension) - 1;
452       index(dimension - 1) = index(dimension - 1) - 1;
453       --dimension;
454     }
455   }
456 }
457
458 template <typename TA, typename TB, typename TC>
459 void MultiplyElements(const ArrayND<TA, 3> &a,
460                       const ArrayND<TB, 3> &b,
461                       ArrayND<TC, 3> *c) {
462   // Specialization for N==3
463   c->ResizeLike(a);
464   assert(a.Shape(0) == b.Shape(0));
465   assert(a.Shape(1) == b.Shape(1));
466   assert(a.Shape(2) == b.Shape(2));
467   for (int i = 0; i < a.Shape(0); ++i) {
468     for (int j = 0; j < a.Shape(1); ++j) {
469       for (int k = 0; k < a.Shape(2); ++k) {
470         (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k));
471       }
472     }
473   }
474 }
475
476 template <typename TA, typename TB, typename TC>
477 void MultiplyElements(const Array3D<TA> &a,
478                       const Array3D<TB> &b,
479                       Array3D<TC> *c) {
480   // Specialization for N==3
481   c->ResizeLike(a);
482   assert(a.Shape(0) == b.Shape(0));
483   assert(a.Shape(1) == b.Shape(1));
484   assert(a.Shape(2) == b.Shape(2));
485   for (int i = 0; i < a.Shape(0); ++i) {
486     for (int j = 0; j < a.Shape(1); ++j) {
487       for (int k = 0; k < a.Shape(2); ++k) {
488         (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k));
489       }
490     }
491   }
492 }
493
494 }  // namespace libmv
495
496 #endif  // LIBMV_IMAGE_ARRAY_ND_H