diff options
Diffstat (limited to 'intern/libmv/libmv/image')
| -rw-r--r-- | intern/libmv/libmv/image/array_nd.cc | 108 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/array_nd.h | 497 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/array_nd_test.cc | 324 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/convolve.cc | 344 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/convolve.h | 107 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/convolve_test.cc | 110 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/correlation.h | 74 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/image.h | 155 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/image_converter.h | 81 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/image_drawing.h | 285 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/image_test.cc | 45 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/sample.h | 138 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/sample_test.cc | 89 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/tuple.h | 90 | ||||
| -rw-r--r-- | intern/libmv/libmv/image/tuple_test.cc | 83 |
15 files changed, 2530 insertions, 0 deletions
diff --git a/intern/libmv/libmv/image/array_nd.cc b/intern/libmv/libmv/image/array_nd.cc new file mode 100644 index 00000000000..469a19aabf1 --- /dev/null +++ b/intern/libmv/libmv/image/array_nd.cc @@ -0,0 +1,108 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include "libmv/image/image.h" +#include <iostream> +#include <cmath> + +namespace libmv { + +void FloatArrayToScaledByteArray(const Array3Df &float_array, + Array3Du *byte_array, + bool automatic_range_detection + ) { + byte_array->ResizeLike(float_array); + float minval = HUGE_VAL; + float maxval = -HUGE_VAL; + if (automatic_range_detection) { + for (int i = 0; i < float_array.Height(); ++i) { + for (int j = 0; j < float_array.Width(); ++j) { + for (int k = 0; k < float_array.Depth(); ++k) { + minval = std::min(minval, float_array(i, j, k)); + maxval = std::max(maxval, float_array(i, j, k)); + } + } + } + } else { + minval = 0; + maxval = 1; + } + for (int i = 0; i < float_array.Height(); ++i) { + for (int j = 0; j < float_array.Width(); ++j) { + for (int k = 0; k < float_array.Depth(); ++k) { + float unscaled = (float_array(i, j, k) - minval) / (maxval - minval); + (*byte_array)(i, j, k) = (unsigned char)(255 * unscaled); + } + } + } +} + +void ByteArrayToScaledFloatArray(const Array3Du &byte_array, + Array3Df *float_array) { + float_array->ResizeLike(byte_array); + for (int i = 0; i < byte_array.Height(); ++i) { + for (int j = 0; j < byte_array.Width(); ++j) { + for (int k = 0; k < byte_array.Depth(); ++k) { + (*float_array)(i, j, k) = float(byte_array(i, j, k)) / 255.0f; + } + } + } +} + +void SplitChannels(const Array3Df &input, + Array3Df *channel0, + Array3Df *channel1, + Array3Df *channel2) { + assert(input.Depth() >= 3); + channel0->Resize(input.Height(), input.Width()); + channel1->Resize(input.Height(), input.Width()); + channel2->Resize(input.Height(), input.Width()); + for (int row = 0; row < input.Height(); ++row) { + for (int column = 0; column < input.Width(); ++column) { + (*channel0)(row, column) = input(row, column, 0); + (*channel1)(row, column) = input(row, column, 1); + (*channel2)(row, column) = input(row, column, 2); + } + } +} + +void PrintArray(const Array3Df &array) { + using namespace std; + + printf("[\n"); + for (int r = 0; r < array.Height(); ++r) { + printf("["); + for (int c = 0; c < array.Width(); ++c) { + if (array.Depth() == 1) { + printf("%11f, ", array(r, c)); + } else { + printf("["); + for (int k = 0; k < array.Depth(); ++k) { + printf("%11f, ", array(r, c, k)); + } + printf("],"); + } + } + printf("],\n"); + } + printf("]\n"); +} + +} // namespace libmv diff --git a/intern/libmv/libmv/image/array_nd.h b/intern/libmv/libmv/image/array_nd.h new file mode 100644 index 00000000000..e95e66aa2b3 --- /dev/null +++ b/intern/libmv/libmv/image/array_nd.h @@ -0,0 +1,497 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_ARRAY_ND_H +#define LIBMV_IMAGE_ARRAY_ND_H + +#include <cassert> +#include <cstdio> +#include <cstring> + +#include "libmv/image/tuple.h" + +namespace libmv { + +class BaseArray {}; + +/// A multidimensional array class. +template <typename T, int N> +class ArrayND : public BaseArray { + public: + typedef T Scalar; + + /// Type for the multidimensional indices. + typedef Tuple<int, N> Index; + + /// Create an empty array. + ArrayND() : data_(NULL), own_data_(true) { Resize(Index(0)); } + + /// Create an array with the specified shape. + ArrayND(const Index &shape) : data_(NULL), own_data_(true) { Resize(shape); } + + /// Create an array with the specified shape. + ArrayND(int *shape) : data_(NULL), own_data_(true) { Resize(shape); } + + /// Copy constructor. + ArrayND(const ArrayND<T, N> &b) : data_(NULL), own_data_(true) { + ResizeLike(b); + std::memcpy(Data(), b.Data(), sizeof(T) * Size()); + } + + ArrayND(int s0) : data_(NULL), own_data_(true) { Resize(s0); } + ArrayND(int s0, int s1) : data_(NULL), own_data_(true) { Resize(s0, s1); } + ArrayND(int s0, int s1, int s2) : data_(NULL), own_data_(true) { + Resize(s0, s1, s2); + } + + ArrayND(T* data, int s0, int s1, int s2) + : shape_(0), strides_(0), data_(data), own_data_(false) { + Resize(s0, s1, s2); + } + + /// Destructor deletes pixel data. + ~ArrayND() { + if (own_data_) { + delete [] data_; + } + } + + /// Assignation copies pixel data. + ArrayND &operator=(const ArrayND<T, N> &b) { + assert(this != &b); + ResizeLike(b); + std::memcpy(Data(), b.Data(), sizeof(T) * Size()); + return *this; + } + + const Index &Shapes() const { + return shape_; + } + + const Index &Strides() const { + return strides_; + } + + /// Create an array of shape s. + void Resize(const Index &new_shape) { + if (data_ != NULL && shape_ == new_shape) { + // Don't bother realloacting if the shapes match. + return; + } + shape_.Reset(new_shape); + strides_(N - 1) = 1; + for (int i = N - 1; i > 0; --i) { + strides_(i - 1) = strides_(i) * shape_(i); + } + if (own_data_) { + delete [] data_; + data_ = NULL; + if (Size() > 0) { + data_ = new T[Size()]; + } + } + } + + template<typename D> + void ResizeLike(const ArrayND<D, N> &other) { + Resize(other.Shape()); + } + + /// Resizes the array to shape s. All data is lost. + void Resize(const int *new_shape_array) { + Resize(Index(new_shape_array)); + } + + /// Resize a 1D array to length s0. + void Resize(int s0) { + assert(N == 1); + int shape[] = {s0}; + Resize(shape); + } + + /// Resize a 2D array to shape (s0,s1). + void Resize(int s0, int s1) { + int shape[N] = {s0, s1}; + for (int i = 2; i < N; ++i) { + shape[i] = 1; + } + Resize(shape); + } + + // Match Eigen2's API. + void resize(int rows, int cols) { + Resize(rows, cols); + } + + /// Resize a 3D array to shape (s0,s1,s2). + void Resize(int s0, int s1, int s2) { + assert(N == 3); + int shape[] = {s0, s1, s2}; + Resize(shape); + } + + template<typename D> + void CopyFrom(const ArrayND<D, N> &other) { + ResizeLike(other); + T *data = Data(); + const D *other_data = other.Data(); + for (int i = 0; i < Size(); ++i) { + data[i] = T(other_data[i]); + } + } + + void Fill(T value) { + for (int i = 0; i < Size(); ++i) { + Data()[i] = value; + } + } + + // Match Eigen's API. + void fill(T value) { + for (int i = 0; i < Size(); ++i) { + Data()[i] = value; + } + } + + /// Return a tuple containing the length of each axis. + const Index &Shape() const { + return shape_; + } + + /// Return the length of an axis. + int Shape(int axis) const { + return shape_(axis); + } + + /// Return the distance between neighboring elements along axis. + int Stride(int axis) const { + return strides_(axis); + } + + /// Return the number of elements of the array. + int Size() const { + int size = 1; + for (int i = 0; i < N; ++i) + size *= Shape(i); + return size; + } + + /// Return the total amount of memory used by the array. + int MemorySizeInBytes() const { + return sizeof(*this) + Size() * sizeof(T); + } + + /// Pointer to the first element of the array. + T *Data() { return data_; } + + /// Constant pointer to the first element of the array. + const T *Data() const { return data_; } + + /// Distance between the first element and the element at position index. + int Offset(const Index &index) const { + int offset = 0; + for (int i = 0; i < N; ++i) + offset += index(i) * Stride(i); + return offset; + } + + /// 1D specialization. + int Offset(int i0) const { + assert(N == 1); + return i0 * Stride(0); + } + + /// 2D specialization. + int Offset(int i0, int i1) const { + assert(N == 2); + return i0 * Stride(0) + i1 * Stride(1); + } + + /// 3D specialization. + int Offset(int i0, int i1, int i2) const { + assert(N == 3); + return i0 * Stride(0) + i1 * Stride(1) + i2 * Stride(2); + } + + /// Return a reference to the element at position index. + T &operator()(const Index &index) { + // TODO(pau) Boundary checking in debug mode. + return *( Data() + Offset(index) ); + } + + /// 1D specialization. + T &operator()(int i0) { + return *( Data() + Offset(i0) ); + } + + /// 2D specialization. + T &operator()(int i0, int i1) { + assert(0 <= i0 && i0 < Shape(0)); + assert(0 <= i1 && i1 < Shape(1)); + return *(Data() + Offset(i0, i1)); + } + + /// 3D specialization. + T &operator()(int i0, int i1, int i2) { + assert(0 <= i0 && i0 < Shape(0)); + assert(0 <= i1 && i1 < Shape(1)); + assert(0 <= i2 && i2 < Shape(2)); + return *(Data() + Offset(i0, i1, i2)); + } + + /// Return a constant reference to the element at position index. + const T &operator()(const Index &index) const { + return *(Data() + Offset(index)); + } + + /// 1D specialization. + const T &operator()(int i0) const { + return *(Data() + Offset(i0)); + } + + /// 2D specialization. + const T &operator()(int i0, int i1) const { + assert(0 <= i0 && i0 < Shape(0)); + assert(0 <= i1 && i1 < Shape(1)); + return *(Data() + Offset(i0, i1)); + } + + /// 3D specialization. + const T &operator()(int i0, int i1, int i2) const { + return *(Data() + Offset(i0, i1, i2)); + } + + /// True if index is inside array. + bool Contains(const Index &index) const { + for (int i = 0; i < N; ++i) + if (index(i) < 0 || index(i) >= Shape(i)) + return false; + return true; + } + + /// 1D specialization. + bool Contains(int i0) const { + return 0 <= i0 && i0 < Shape(0); + } + + /// 2D specialization. + bool Contains(int i0, int i1) const { + return 0 <= i0 && i0 < Shape(0) + && 0 <= i1 && i1 < Shape(1); + } + + /// 3D specialization. + bool Contains(int i0, int i1, int i2) const { + return 0 <= i0 && i0 < Shape(0) + && 0 <= i1 && i1 < Shape(1) + && 0 <= i2 && i2 < Shape(2); + } + + bool operator==(const ArrayND<T, N> &other) const { + if (shape_ != other.shape_) return false; + if (strides_ != other.strides_) return false; + for (int i = 0; i < Size(); ++i) { + if (this->Data()[i] != other.Data()[i]) + return false; + } + return true; + } + + bool operator!=(const ArrayND<T, N> &other) const { + return !(*this == other); + } + + ArrayND<T, N> operator*(const ArrayND<T, N> &other) const { + assert(Shape() = other.Shape()); + ArrayND<T, N> res; + res.ResizeLike(*this); + for (int i = 0; i < res.Size(); ++i) { + res.Data()[i] = Data()[i] * other.Data()[i]; + } + return res; + } + + protected: + /// The number of element in each dimension. + Index shape_; + + /// How to jump to neighbors in each dimension. + Index strides_; + + /// Pointer to the first element of the array. + T *data_; + + /// Flag if this Array either own or reference the data + bool own_data_; +}; + +/// 3D array (row, column, channel). +template <typename T> +class Array3D : public ArrayND<T, 3> { + typedef ArrayND<T, 3> Base; + public: + Array3D() + : Base() { + } + Array3D(int height, int width, int depth = 1) + : Base(height, width, depth) { + } + Array3D(T* data, int height, int width, int depth = 1) + : Base(data, height, width, depth) { + } + + void Resize(int height, int width, int depth = 1) { + Base::Resize(height, width, depth); + } + + int Height() const { + return Base::Shape(0); + } + int Width() const { + return Base::Shape(1); + } + int Depth() const { + return Base::Shape(2); + } + + // Match Eigen2's API so that Array3D's and Mat*'s can work together via + // template magic. + int rows() const { return Height(); } + int cols() const { return Width(); } + int depth() const { return Depth(); } + + int Get_Step() const { return Width()*Depth(); } + + /// Enable accessing with 2 indices for grayscale images. + T &operator()(int i0, int i1, int i2 = 0) { + assert(0 <= i0 && i0 < Height()); + assert(0 <= i1 && i1 < Width()); + return Base::operator()(i0, i1, i2); + } + const T &operator()(int i0, int i1, int i2 = 0) const { + assert(0 <= i0 && i0 < Height()); + assert(0 <= i1 && i1 < Width()); + return Base::operator()(i0, i1, i2); + } +}; + +typedef Array3D<unsigned char> Array3Du; +typedef Array3D<unsigned int> Array3Dui; +typedef Array3D<int> Array3Di; +typedef Array3D<float> Array3Df; +typedef Array3D<short> Array3Ds; + +void SplitChannels(const Array3Df &input, + Array3Df *channel0, + Array3Df *channel1, + Array3Df *channel2); + +void PrintArray(const Array3Df &array); + +/** Convert a float array into a byte array by scaling values by 255* (max-min). + * where max and min are automatically detected + * (if automatic_range_detection = true) + * \note and TODO this automatic detection only works when the image contains + * at least one pixel of both bounds. + **/ +void FloatArrayToScaledByteArray(const Array3Df &float_array, + Array3Du *byte_array, + bool automatic_range_detection = false); + +//! Convert a byte array into a float array by dividing values by 255. +void ByteArrayToScaledFloatArray(const Array3Du &byte_array, + Array3Df *float_array); + +template <typename AArrayType, typename BArrayType, typename CArrayType> +void MultiplyElements(const AArrayType &a, + const BArrayType &b, + CArrayType *c) { + // This function does an element-wise multiply between + // the two Arrays A and B, and stores the result in C. + // A and B must have the same dimensions. + assert(a.Shape() == b.Shape()); + c->ResizeLike(a); + + // To perform the multiplcation, a "current" index into the N-dimensions of + // the A and B matrix specifies which elements are being multiplied. + typename CArrayType::Index index; + + // The index starts at the maximum value for each dimension + const typename CArrayType::Index& cShape = c->Shape(); + for ( int i = 0; i < CArrayType::Index::SIZE; ++i ) + index(i) = cShape(i) - 1; + + // After each multiplication, the highest-dimensional index is reduced. + // if this reduces it less than zero, it resets to its maximum value + // and decrements the index of the next lower dimension. + // This ripple-action continues until the entire new array has been + // calculated, indicated by dimension zero having a negative index. + while ( index(0) >= 0 ) { + (*c)(index) = a(index) * b(index); + + int dimension = CArrayType::Index::SIZE - 1; + index(dimension) = index(dimension) - 1; + while ( dimension > 0 && index(dimension) < 0 ) { + index(dimension) = cShape(dimension) - 1; + index(dimension - 1) = index(dimension - 1) - 1; + --dimension; + } + } +} + +template <typename TA, typename TB, typename TC> +void MultiplyElements(const ArrayND<TA, 3> &a, + const ArrayND<TB, 3> &b, + ArrayND<TC, 3> *c) { + // Specialization for N==3 + c->ResizeLike(a); + assert(a.Shape(0) == b.Shape(0)); + assert(a.Shape(1) == b.Shape(1)); + assert(a.Shape(2) == b.Shape(2)); + for (int i = 0; i < a.Shape(0); ++i) { + for (int j = 0; j < a.Shape(1); ++j) { + for (int k = 0; k < a.Shape(2); ++k) { + (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k)); + } + } + } +} + +template <typename TA, typename TB, typename TC> +void MultiplyElements(const Array3D<TA> &a, + const Array3D<TB> &b, + Array3D<TC> *c) { + // Specialization for N==3 + c->ResizeLike(a); + assert(a.Shape(0) == b.Shape(0)); + assert(a.Shape(1) == b.Shape(1)); + assert(a.Shape(2) == b.Shape(2)); + for (int i = 0; i < a.Shape(0); ++i) { + for (int j = 0; j < a.Shape(1); ++j) { + for (int k = 0; k < a.Shape(2); ++k) { + (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k)); + } + } + } +} + +} // namespace libmv + +#endif // LIBMV_IMAGE_ARRAY_ND_H diff --git a/intern/libmv/libmv/image/array_nd_test.cc b/intern/libmv/libmv/image/array_nd_test.cc new file mode 100644 index 00000000000..313f21b60e9 --- /dev/null +++ b/intern/libmv/libmv/image/array_nd_test.cc @@ -0,0 +1,324 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include "libmv/image/array_nd.h" +#include "testing/testing.h" + +using libmv::ArrayND; +using libmv::Array3D; +using libmv::Array3Df; + +namespace { + +TEST(ArrayND, EmptyConstructor) { + ArrayND<int, 2> a; + + EXPECT_EQ(0, a.Shape(0)); + EXPECT_EQ(0, a.Shape(1)); +} + +TEST(ArrayND, IndexConstructor) { + int s[] = {1, 2, 3}; + ArrayND<int, 3>::Index shape(s); + ArrayND<int, 3> a(shape); + + EXPECT_EQ(1, a.Shape(0)); + EXPECT_EQ(2, a.Shape(1)); + EXPECT_EQ(3, a.Shape(2)); +} + +TEST(ArrayND, PointerConstructor) { + int s[] = {1, 2, 3}; + ArrayND<int, 3> a(s); + + EXPECT_EQ(1, a.Shape(0)); + EXPECT_EQ(2, a.Shape(1)); + EXPECT_EQ(3, a.Shape(2)); +} + +TEST(ArrayND, CopyConstructor) { + int s[] = {1, 2, 3}; + ArrayND<int, 3> a(s); + a(0, 1, 1) = 3; + a(0, 1, 2) = 3; + ArrayND<int, 3> b(a); + EXPECT_EQ(1, b.Shape(0)); + EXPECT_EQ(2, b.Shape(1)); + EXPECT_EQ(3, b.Shape(2)); + EXPECT_EQ(3, b(0, 1, 1)); + b(0, 1, 2) = 2; + EXPECT_EQ(3, a(0, 1, 2)); +} + +TEST(ArrayND, Assignation) { + int s[] = {1, 2, 3}; + ArrayND<int, 3> a(s); + a(0, 1, 1) = 3; + a(0, 1, 2) = 3; + ArrayND<int, 3> b; + b = a; + EXPECT_EQ(1, b.Shape(0)); + EXPECT_EQ(2, b.Shape(1)); + EXPECT_EQ(3, b.Shape(2)); + EXPECT_EQ(3, b(0, 1, 1)); + b(0, 1, 2) = 2; + EXPECT_EQ(3, a(0, 1, 2)); +} + +TEST(ArrayND, Fill) { + int s[] = {2, 2}; + ArrayND<int, 2> a(s); + a.Fill(42); + EXPECT_EQ(42, a(0, 0)); + EXPECT_EQ(42, a(0, 1)); + EXPECT_EQ(42, a(1, 0)); + EXPECT_EQ(42, a(1, 1)); +} + +TEST(ArrayND, Size) { + int s[] = {1, 2, 3}; + ArrayND<int, 3>::Index shape(s); + ArrayND<int, 3> a(shape); + + int l[] = {0, 1, 2}; + ArrayND<int, 3>::Index last(l); + + EXPECT_EQ(a.Size(), a.Offset(last)+1); + EXPECT_TRUE(a.Contains(last)); + EXPECT_FALSE(a.Contains(shape)); +} + +TEST(ArrayND, MemorySizeInBytes) { + int s[] = {2, 3}; + ArrayND<int, 2>::Index shape(s); + ArrayND<int, 2> a(shape); + + int size = 24 + sizeof(a); + EXPECT_EQ(size, a.MemorySizeInBytes()); +} + +TEST(ArrayND, Parenthesis) { + typedef ArrayND<int, 2>::Index Index; + + int s[] = {3, 3}; + ArrayND<int, 2> a(s); + + *(a.Data()+0) = 0; + *(a.Data()+5) = 5; + + int i1[] = {0, 0}; + EXPECT_EQ(0, a(Index(i1))); + int i2[] = {1, 2}; + EXPECT_EQ(5, a(Index(i2))); +} + +TEST(ArrayND, 1DConstructor) { + ArrayND<int, 1> a(3); + + EXPECT_EQ(3, a.Shape(0)); +} + +TEST(ArrayND, 2DConstructor) { + ArrayND<int, 2> a(1, 2); + + EXPECT_EQ(1, a.Shape(0)); + EXPECT_EQ(2, a.Shape(1)); +} + +TEST(ArrayND, 3DConstructor) { + ArrayND<int, 3> a(1, 2, 3); + + EXPECT_EQ(1, a.Shape(0)); + EXPECT_EQ(2, a.Shape(1)); + EXPECT_EQ(3, a.Shape(2)); +} + +TEST(ArrayND, 1DAccessor) { + ArrayND<int, 1> a(3); + a(0) = 1; + a(1) = 2; + + EXPECT_EQ(1, a(0)); + EXPECT_EQ(2, a(1)); + EXPECT_EQ(1, *(a.Data())); + EXPECT_EQ(2, *(a.Data() + a.Stride(0))); +} + +TEST(ArrayND, 2DAccessor) { + ArrayND<int, 2> a(3, 3); + a(0, 0) = 1; + a(1, 1) = 2; + + EXPECT_EQ(1, a(0, 0)); + EXPECT_EQ(2, a(1, 1)); + EXPECT_EQ(1, *(a.Data())); + EXPECT_EQ(2, *(a.Data() + a.Stride(0) + a.Stride(1))); +} + +TEST(ArrayND, 3DAccessor) { + ArrayND<int, 3> a(3, 3, 3); + a(0, 0, 0) = 1; + a(1, 1, 1) = 2; + + EXPECT_EQ(1, a(0, 0, 0)); + EXPECT_EQ(2, a(1, 1, 1)); + EXPECT_EQ(1, *(a.Data())); + EXPECT_EQ(2, *(a.Data() + a.Stride(0) + a.Stride(1) + a.Stride(2))); +} + +TEST(ArrayND, CopyFrom) { + ArrayND<int, 3> a(2, 2, 1); + a(0, 0, 0) = 1; + a(0, 1, 0) = 2; + a(1, 0, 0) = 3; + a(1, 1, 0) = 4; + ArrayND<float, 3> b; + b.CopyFrom(a); + EXPECT_FLOAT_EQ(1.f, b(0, 0, 0)); + EXPECT_FLOAT_EQ(2.f, b(0, 1, 0)); + EXPECT_FLOAT_EQ(3.f, b(1, 0, 0)); + EXPECT_FLOAT_EQ(4.f, b(1, 1, 0)); +} + +TEST(ArrayND, MultiplyElements) { + ArrayND<int, 3> a(2, 2, 1); + a(0, 0, 0) = 1; + a(0, 1, 0) = 2; + a(1, 0, 0) = 3; + a(1, 1, 0) = 4; + ArrayND<int, 3> b(2, 2, 1); + b(0, 0, 0) = 6; + b(0, 1, 0) = 5; + b(1, 0, 0) = 4; + b(1, 1, 0) = 3; + ArrayND<int, 3> c; + MultiplyElements(a, b, &c); + EXPECT_FLOAT_EQ(6, c(0, 0, 0)); + EXPECT_FLOAT_EQ(10, c(0, 1, 0)); + EXPECT_FLOAT_EQ(12, c(1, 0, 0)); + EXPECT_FLOAT_EQ(12, c(1, 1, 0)); +} + +TEST(ArrayND, IsEqualOperator) { + ArrayND<int, 3> a(2, 2, 1); + a(0, 0, 0) = 1; + a(0, 1, 0) = 2; + a(1, 0, 0) = 3; + a(1, 1, 0) = 4; + ArrayND<int, 3> b(2, 2, 1); + b(0, 0, 0) = 1; + b(0, 1, 0) = 2; + b(1, 0, 0) = 3; + b(1, 1, 0) = 4; + EXPECT_TRUE(a == b); + EXPECT_FALSE(a != b); + b(1, 1, 0) = 5; + EXPECT_TRUE(a != b); + EXPECT_FALSE(a == b); +} + +TEST(Array3D, Sizes) { + Array3D<int> array; + EXPECT_EQ(array.Height(), 0); + EXPECT_EQ(array.Width(), 0); + EXPECT_EQ(array.Depth(), 0); + EXPECT_EQ(array.Shape(0), 0); +} + +TEST(Array3D, CopyConstructor) { + Array3D<int> array(10, 10); + array(0, 0) = 1; + array(0, 1) = 1; + Array3D<int> copy(array); + EXPECT_EQ(copy.Height(), 10); + EXPECT_EQ(copy.Width(), 10); + EXPECT_EQ(copy.Depth(), 1); + EXPECT_EQ(copy(0, 0), 1); + copy(0, 1) = 2; + EXPECT_EQ(array(0, 1), 1); +} + +TEST(Array3D, Assignation) { + Array3D<int> array(10, 10); + array(0, 0) = 1; + array(0, 1) = 1; + Array3D<int> copy; + copy = array; + EXPECT_EQ(copy.Height(), 10); + EXPECT_EQ(copy.Width(), 10); + EXPECT_EQ(copy.Depth(), 1); + EXPECT_EQ(copy(0, 0), 1); + copy(0, 1) = 2; + EXPECT_EQ(array(0, 1), 1); +} + +TEST(Array3D, Parenthesis) { + Array3D<int> array(1, 2, 3); + array(0, 1, 0) = 3; + EXPECT_EQ(array(0, 1), 3); +} + +TEST(Array3Df, SplitChannels) { + Array3Df array(1, 2, 3); + array(0, 0, 0) = 1; + array(0, 1, 0) = 1; + array(0, 0, 1) = 2; + array(0, 1, 1) = 2; + array(0, 0, 2) = 3; + array(0, 1, 2) = 3; + Array3Df c0, c1, c2; + SplitChannels(array, &c0, &c1, &c2); + for (int row = 0; row < 1; ++row) { + for (int column = 0; column < 2; ++column) { + EXPECT_EQ(array(row, column, 0), c0(row, column)); + EXPECT_EQ(array(row, column, 1), c1(row, column)); + EXPECT_EQ(array(row, column, 2), c2(row, column)); + } + } +} + +TEST(ArrayND, MultiplyElementsGeneric) { + ArrayND<double, 5> A; + ArrayND<int, 5> B; + ArrayND<double, 5> C; + int shape[] = {1, 3, 5, 7, 1}; + A.Resize(shape); + B.Resize(shape); + + A.Fill(1.1); + B.Fill(2); + MultiplyElements(A, B, &C); + + ArrayND<double, 5>::Index cIndex; + for (int d0 = 0; d0 < shape[0]; ++d0) + for (int d1 = 0; d1 < shape[1]; ++d1) + for (int d2 = 0; d2 < shape[2]; ++d2) + for (int d3 = 0; d3 < shape[3]; ++d3) + for (int d4 = 0; d4 < shape[4]; ++d4) { + cIndex(0) = d0; + cIndex(1) = d1; + cIndex(2) = d2; + cIndex(3) = d3; + cIndex(4) = d4; + EXPECT_EQ(2.2, C(cIndex)); + } +} + +} // namespace diff --git a/intern/libmv/libmv/image/convolve.cc b/intern/libmv/libmv/image/convolve.cc new file mode 100644 index 00000000000..464043581d2 --- /dev/null +++ b/intern/libmv/libmv/image/convolve.cc @@ -0,0 +1,344 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include "libmv/image/convolve.h" + +#include <cmath> + +#include "libmv/image/image.h" + +namespace libmv { + +// Compute a Gaussian kernel and derivative, such that you can take the +// derivative of an image by convolving with the kernel horizontally then the +// derivative vertically to get (eg) the y derivative. +void ComputeGaussianKernel(double sigma, Vec *kernel, Vec *derivative) { + assert(sigma >= 0.0); + + // 0.004 implies a 3 pixel kernel with 1 pixel sigma. + const float truncation_factor = 0.004f; + + // Calculate the kernel size based on sigma such that it is odd. + float precisehalfwidth = GaussianInversePositive(truncation_factor, sigma); + int width = lround(2*precisehalfwidth); + if (width % 2 == 0) { + width++; + } + // Calculate the gaussian kernel and its derivative. + kernel->resize(width); + derivative->resize(width); + kernel->setZero(); + derivative->setZero(); + int halfwidth = width / 2; + for (int i = -halfwidth; i <= halfwidth; ++i) { + (*kernel)(i + halfwidth) = Gaussian(i, sigma); + (*derivative)(i + halfwidth) = GaussianDerivative(i, sigma); + } + // Since images should not get brighter or darker, normalize. + NormalizeL1(kernel); + + // Normalize the derivative differently. See + // www.cs.duke.edu/courses/spring03/cps296.1/handouts/Image%20Processing.pdf + double factor = 0.; + for (int i = -halfwidth; i <= halfwidth; ++i) { + factor -= i*(*derivative)(i+halfwidth); + } + *derivative /= factor; +} + +template <int size, bool vertical> +void FastConvolve(const Vec &kernel, int width, int height, + const float* src, int src_stride, int src_line_stride, + float* dst, int dst_stride) { + double coefficients[2 * size + 1]; + for (int k = 0; k < 2 * size + 1; ++k) { + coefficients[k] = kernel(2 * size - k); + } + // Fast path: if the kernel has a certain size, use the constant sized loops. + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + double sum = 0; + for (int k = -size; k <= size; ++k) { + if (vertical) { + if (y + k >= 0 && y + k < height) { + sum += src[k * src_line_stride] * coefficients[k + size]; + } + } else { + if (x + k >= 0 && x + k < width) { + sum += src[k * src_stride] * coefficients[k + size]; + } + } + } + dst[0] = static_cast<float>(sum); + src += src_stride; + dst += dst_stride; + } + } +} + +template<bool vertical> +void Convolve(const Array3Df &in, + const Vec &kernel, + Array3Df *out_pointer, + int plane) { + int width = in.Width(); + int height = in.Height(); + Array3Df &out = *out_pointer; + if (plane == -1) { + out.ResizeLike(in); + plane = 0; + } + + assert(kernel.size() % 2 == 1); + assert(&in != out_pointer); + + int src_line_stride = in.Stride(0); + int src_stride = in.Stride(1); + int dst_stride = out.Stride(1); + const float* src = in.Data(); + float* dst = out.Data() + plane; + + // Use a dispatch table to make most convolutions used in practice use the + // fast path. + int half_width = kernel.size() / 2; + switch (half_width) { +#define static_convolution(size) case size: \ + FastConvolve<size, vertical>(kernel, width, height, src, src_stride, \ + src_line_stride, dst, dst_stride); break; + static_convolution(1) + static_convolution(2) + static_convolution(3) + static_convolution(4) + static_convolution(5) + static_convolution(6) + static_convolution(7) +#undef static_convolution + default: + int dynamic_size = kernel.size() / 2; + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + double sum = 0; + // Slow path: this loop cannot be unrolled. + for (int k = -dynamic_size; k <= dynamic_size; ++k) { + if (vertical) { + if (y + k >= 0 && y + k < height) { + sum += src[k * src_line_stride] * + kernel(2 * dynamic_size - (k + dynamic_size)); + } + } else { + if (x + k >= 0 && x + k < width) { + sum += src[k * src_stride] * + kernel(2 * dynamic_size - (k + dynamic_size)); + } + } + } + dst[0] = static_cast<float>(sum); + src += src_stride; + dst += dst_stride; + } + } + } +} + +void ConvolveHorizontal(const Array3Df &in, + const Vec &kernel, + Array3Df *out_pointer, + int plane) { + Convolve<false>(in, kernel, out_pointer, plane); +} + +void ConvolveVertical(const Array3Df &in, + const Vec &kernel, + Array3Df *out_pointer, + int plane) { + Convolve<true>(in, kernel, out_pointer, plane); +} + +void ConvolveGaussian(const Array3Df &in, + double sigma, + Array3Df *out_pointer) { + Vec kernel, derivative; + ComputeGaussianKernel(sigma, &kernel, &derivative); + + Array3Df tmp; + ConvolveVertical(in, kernel, &tmp); + ConvolveHorizontal(tmp, kernel, out_pointer); +} + +void ImageDerivatives(const Array3Df &in, + double sigma, + Array3Df *gradient_x, + Array3Df *gradient_y) { + Vec kernel, derivative; + ComputeGaussianKernel(sigma, &kernel, &derivative); + Array3Df tmp; + + // Compute first derivative in x. + ConvolveVertical(in, kernel, &tmp); + ConvolveHorizontal(tmp, derivative, gradient_x); + + // Compute first derivative in y. + ConvolveHorizontal(in, kernel, &tmp); + ConvolveVertical(tmp, derivative, gradient_y); +} + +void BlurredImageAndDerivatives(const Array3Df &in, + double sigma, + Array3Df *blurred_image, + Array3Df *gradient_x, + Array3Df *gradient_y) { + Vec kernel, derivative; + ComputeGaussianKernel(sigma, &kernel, &derivative); + Array3Df tmp; + + // Compute convolved image. + ConvolveVertical(in, kernel, &tmp); + ConvolveHorizontal(tmp, kernel, blurred_image); + + // Compute first derivative in x (reusing vertical convolution above). + ConvolveHorizontal(tmp, derivative, gradient_x); + + // Compute first derivative in y. + ConvolveHorizontal(in, kernel, &tmp); + ConvolveVertical(tmp, derivative, gradient_y); +} + +// Compute the gaussian blur of an image and the derivatives of the blurred +// image, and store the results in three channels. Since the blurred value and +// gradients are closer in memory, this leads to better performance if all +// three values are needed at the same time. +void BlurredImageAndDerivativesChannels(const Array3Df &in, + double sigma, + Array3Df *blurred_and_gradxy) { + assert(in.Depth() == 1); + + Vec kernel, derivative; + ComputeGaussianKernel(sigma, &kernel, &derivative); + + // Compute convolved image. + Array3Df tmp; + ConvolveVertical(in, kernel, &tmp); + blurred_and_gradxy->Resize(in.Height(), in.Width(), 3); + ConvolveHorizontal(tmp, kernel, blurred_and_gradxy, 0); + + // Compute first derivative in x. + ConvolveHorizontal(tmp, derivative, blurred_and_gradxy, 1); + + // Compute first derivative in y. + ConvolveHorizontal(in, kernel, &tmp); + ConvolveVertical(tmp, derivative, blurred_and_gradxy, 2); +} + +void BoxFilterHorizontal(const Array3Df &in, + int window_size, + Array3Df *out_pointer) { + Array3Df &out = *out_pointer; + out.ResizeLike(in); + int half_width = (window_size - 1) / 2; + + for (int k = 0; k < in.Depth(); ++k) { + for (int i = 0; i < in.Height(); ++i) { + float sum = 0; + // Init sum. + for (int j = 0; j < half_width; ++j) { + sum += in(i, j, k); + } + // Fill left border. + for (int j = 0; j < half_width + 1; ++j) { + sum += in(i, j + half_width, k); + out(i, j, k) = sum; + } + // Fill interior. + for (int j = half_width + 1; j < in.Width()-half_width; ++j) { + sum -= in(i, j - half_width - 1, k); + sum += in(i, j + half_width, k); + out(i, j, k) = sum; + } + // Fill right border. + for (int j = in.Width() - half_width; j < in.Width(); ++j) { + sum -= in(i, j - half_width - 1, k); + out(i, j, k) = sum; + } + } + } +} + +void BoxFilterVertical(const Array3Df &in, + int window_size, + Array3Df *out_pointer) { + Array3Df &out = *out_pointer; + out.ResizeLike(in); + int half_width = (window_size - 1) / 2; + + for (int k = 0; k < in.Depth(); ++k) { + for (int j = 0; j < in.Width(); ++j) { + float sum = 0; + // Init sum. + for (int i = 0; i < half_width; ++i) { + sum += in(i, j, k); + } + // Fill left border. + for (int i = 0; i < half_width + 1; ++i) { + sum += in(i + half_width, j, k); + out(i, j, k) = sum; + } + // Fill interior. + for (int i = half_width + 1; i < in.Height()-half_width; ++i) { + sum -= in(i - half_width - 1, j, k); + sum += in(i + half_width, j, k); + out(i, j, k) = sum; + } + // Fill right border. + for (int i = in.Height() - half_width; i < in.Height(); ++i) { + sum -= in(i - half_width - 1, j, k); + out(i, j, k) = sum; + } + } + } +} + +void BoxFilter(const Array3Df &in, + int box_width, + Array3Df *out) { + Array3Df tmp; + BoxFilterHorizontal(in, box_width, &tmp); + BoxFilterVertical(tmp, box_width, out); +} + +void LaplaceFilter(unsigned char* src, + unsigned char* dst, + int width, + int height, + int strength) { + for (int y = 1; y < height-1; y++) + for (int x = 1; x < width-1; x++) { + const unsigned char* s = &src[y*width+x]; + int l = 128 + + s[-width-1] + s[-width] + s[-width+1] + + s[1] - 8*s[0] + s[1] + + s[ width-1] + s[ width] + s[ width+1]; + int d = ((256-strength)*s[0] + strength*l) / 256; + if (d < 0) d=0; + if (d > 255) d=255; + dst[y*width+x] = d; + } +} + +} // namespace libmv diff --git a/intern/libmv/libmv/image/convolve.h b/intern/libmv/libmv/image/convolve.h new file mode 100644 index 00000000000..d3b6da9794b --- /dev/null +++ b/intern/libmv/libmv/image/convolve.h @@ -0,0 +1,107 @@ +// Copyright (c) 2007, 2008, 2011 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_CONVOLVE_H_ +#define LIBMV_IMAGE_CONVOLVE_H_ + +#include "libmv/image/image.h" +#include "libmv/numeric/numeric.h" + +namespace libmv { + +// TODO(keir): Find a better place for these functions. gaussian.h in numeric? + +// Zero mean Gaussian. +inline double Gaussian(double x, double sigma) { + return 1/sqrt(2*M_PI*sigma*sigma) * exp(-(x*x/2/sigma/sigma)); +} +// 2D gaussian (zero mean) +// (9) in http://mathworld.wolfram.com/GaussianFunction.html +inline double Gaussian2D(double x, double y, double sigma) { + return 1.0/(2.0*M_PI*sigma*sigma) * exp( -(x*x+y*y)/(2.0*sigma*sigma)); +} +inline double GaussianDerivative(double x, double sigma) { + return -x / sigma / sigma * Gaussian(x, sigma); +} +// Solve the inverse of the Gaussian for positive x. +inline double GaussianInversePositive(double y, double sigma) { + return sqrt(-2 * sigma * sigma * log(y * sigma * sqrt(2*M_PI))); +} + +void ComputeGaussianKernel(double sigma, Vec *kernel, Vec *derivative); +void ConvolveHorizontal(const FloatImage &in, + const Vec &kernel, + FloatImage *out_pointer, + int plane = -1); +void ConvolveVertical(const FloatImage &in, + const Vec &kernel, + FloatImage *out_pointer, + int plane = -1); +void ConvolveGaussian(const FloatImage &in, + double sigma, + FloatImage *out_pointer); + +void ImageDerivatives(const FloatImage &in, + double sigma, + FloatImage *gradient_x, + FloatImage *gradient_y); + +void BlurredImageAndDerivatives(const FloatImage &in, + double sigma, + FloatImage *blurred_image, + FloatImage *gradient_x, + FloatImage *gradient_y); + +// Blur and take the gradients of an image, storing the results inside the +// three channels of blurred_and_gradxy. +void BlurredImageAndDerivativesChannels(const FloatImage &in, + double sigma, + FloatImage *blurred_and_gradxy); + +void BoxFilterHorizontal(const FloatImage &in, + int window_size, + FloatImage *out_pointer); + +void BoxFilterVertical(const FloatImage &in, + int window_size, + FloatImage *out_pointer); + +void BoxFilter(const FloatImage &in, + int box_width, + FloatImage *out); + +/*! + Convolve \a src into \a dst with the discrete laplacian operator. + + \a src and \a dst should be \a width x \a height images. + \a strength is an interpolation coefficient (0-256) between original image and the laplacian. + + \note Make sure the search region is filtered with the same strength as the pattern. +*/ +void LaplaceFilter(unsigned char* src, + unsigned char* dst, + int width, + int height, + int strength); + +} // namespace libmv + +#endif // LIBMV_IMAGE_CONVOLVE_H_ + diff --git a/intern/libmv/libmv/image/convolve_test.cc b/intern/libmv/libmv/image/convolve_test.cc new file mode 100644 index 00000000000..0cdef8e1e72 --- /dev/null +++ b/intern/libmv/libmv/image/convolve_test.cc @@ -0,0 +1,110 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include <iostream> + +#include "libmv/image/convolve.h" +#include "libmv/image/image.h" +#include "libmv/numeric/numeric.h" +#include "testing/testing.h" + +using namespace libmv; + +namespace { + +TEST(Convolve, ComputeGaussianKernel) { + Vec kernel, derivative; + ComputeGaussianKernel(1, &kernel, &derivative); + EXPECT_EQ(7, kernel.size()); + // TODO(keir): Put in a more thorough test here! +} + +TEST(Convolve, ConvolveGaussian) { + FloatImage im(10, 10); + im.Fill(1); + FloatImage blured; + ConvolveGaussian(im, 3, &blured); + EXPECT_NEAR(im(5, 5), 1, 1e-7); +} + +TEST(Convolve, BoxFilterHorizontal) { + FloatImage im(10, 10), convolved, filtered; + im.Fill(1); + BoxFilterHorizontal(im, 3, &filtered); + Vec kernel(3); + kernel.setConstant(1.); + ConvolveHorizontal(im, kernel, &convolved); + EXPECT_EQ(filtered(5, 5), 3); + EXPECT_TRUE(filtered == convolved); +} + +TEST(Convolve, BoxFilter) { + FloatImage image(5, 5), filtered; + // A single 1.0 inside a 5x5 image should expand to a 3x3 square. + image.Fill(0); + image(2, 2) = 1.0; + BoxFilter(image, 3, &filtered); + for (int j = 0; j < 5; j++) { + for (int i = 0; i < 5; i++) { + if (i == 0 || i == 4 || j == 0 || j == 4) { + EXPECT_EQ(0.0, filtered(j, i)); + } else { + EXPECT_EQ(1.0, filtered(j, i)); + } + } + } +} + +TEST(Convolve, BlurredImageAndDerivativesChannelsFlat) { + FloatImage im(10, 10), blurred_and_derivatives; + im.Fill(1); + BlurredImageAndDerivativesChannels(im, 1.0, &blurred_and_derivatives); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 0), 1.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 1), 0.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 2), 0.0, 1e-7); +} + +TEST(Convolve, BlurredImageAndDerivativesChannelsHorizontalSlope) { + FloatImage image(10, 10), blurred_and_derivatives; + for (int j = 0; j < 10; ++j) { + for (int i = 0; i < 10; ++i) { + image(j, i) = 2*i; + } + } + BlurredImageAndDerivativesChannels(image, 0.9, &blurred_and_derivatives); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 0), 10.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 1), 2.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 2), 0.0, 1e-7); +} + +TEST(Convolve, BlurredImageAndDerivativesChannelsVerticalSlope) { + FloatImage image(10, 10), blurred_and_derivatives; + for (int j = 0; j < 10; ++j) { + for (int i = 0; i < 10; ++i) { + image(j, i) = 2*j; + } + } + BlurredImageAndDerivativesChannels(image, 0.9, &blurred_and_derivatives); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 0), 10.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 1), 0.0, 1e-7); + EXPECT_NEAR(blurred_and_derivatives(5, 5, 2), 2.0, 1e-7); +} + +} // namespace diff --git a/intern/libmv/libmv/image/correlation.h b/intern/libmv/libmv/image/correlation.h new file mode 100644 index 00000000000..c354f7e891e --- /dev/null +++ b/intern/libmv/libmv/image/correlation.h @@ -0,0 +1,74 @@ +// Copyright (c) 2012 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_CORRELATION_H +#define LIBMV_IMAGE_CORRELATION_H + +#include "libmv/logging/logging.h" +#include "libmv/image/image.h" + +namespace libmv { + +inline double PearsonProductMomentCorrelation( + const FloatImage &image_and_gradient1_sampled, + const FloatImage &image_and_gradient2_sampled) { + assert(image_and_gradient1_sampled.Width() == + image_and_gradient2_sampled.Width()); + assert(image_and_gradient1_sampled.Height() == + image_and_gradient2_sampled.Height()); + + const int width = image_and_gradient1_sampled.Width(), + height = image_and_gradient1_sampled.Height(); + double sX = 0, sY = 0, sXX = 0, sYY = 0, sXY = 0; + + for (int r = 0; r < height; ++r) { + for (int c = 0; c < width; ++c) { + double x = image_and_gradient1_sampled(r, c, 0); + double y = image_and_gradient2_sampled(r, c, 0); + sX += x; + sY += y; + sXX += x * x; + sYY += y * y; + sXY += x * y; + } + } + + // Normalize. + double N = width * height; + sX /= N; + sY /= N; + sXX /= N; + sYY /= N; + sXY /= N; + + double var_x = sXX - sX * sX; + double var_y = sYY - sY * sY; + double covariance_xy = sXY - sX * sY; + + double correlation = covariance_xy / sqrt(var_x * var_y); + LG << "Covariance xy: " << covariance_xy + << ", var 1: " << var_x << ", var 2: " << var_y + << ", correlation: " << correlation; + return correlation; +} + +} // namespace libmv + +#endif // LIBMV_IMAGE_IMAGE_CORRELATION_H diff --git a/intern/libmv/libmv/image/image.h b/intern/libmv/libmv/image/image.h new file mode 100644 index 00000000000..e0f200a4c93 --- /dev/null +++ b/intern/libmv/libmv/image/image.h @@ -0,0 +1,155 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_IMAGE_H +#define LIBMV_IMAGE_IMAGE_H + +#include <cmath> + +#include "libmv/image/array_nd.h" + +namespace libmv { + +typedef Array3Du ByteImage; // For backwards compatibility. +typedef Array3Df FloatImage; + +// Type added only to manage special 2D array for feature detection +typedef Array3Di IntImage; +typedef Array3Ds ShortImage; + +// An image class that is a thin wrapper around Array3D's of various types. +// TODO(keir): Decide if we should add reference counting semantics... Maybe it +// is the best solution after all. +class Image { + public: + + // Create an image from an array. The image takes ownership of the array. + Image(Array3Du *array) : array_type_(BYTE), array_(array) {} + Image(Array3Df *array) : array_type_(FLOAT), array_(array) {} + + Image(const Image &img): array_type_(NONE), array_(NULL) { + *this = img; + } + + // Underlying data type. + enum DataType { + NONE, + BYTE, + FLOAT, + INT, + SHORT, + }; + + // Size in bytes that the image takes in memory. + int MemorySizeInBytes() { + int size; + switch (array_type_) { + case BYTE: + size = reinterpret_cast<Array3Du *>(array_)->MemorySizeInBytes(); + break; + case FLOAT: + size = reinterpret_cast<Array3Df *>(array_)->MemorySizeInBytes(); + break; + case INT: + size = reinterpret_cast<Array3Di *>(array_)->MemorySizeInBytes(); + break; + case SHORT: + size = reinterpret_cast<Array3Ds *>(array_)->MemorySizeInBytes(); + break; + default : + size = 0; + assert(0); + } + size += sizeof(*this); + return size; + } + + ~Image() { + switch (array_type_) { + case BYTE: + delete reinterpret_cast<Array3Du *>(array_); + + break; + case FLOAT: + delete reinterpret_cast<Array3Df *>(array_); + + break; + case INT: + delete reinterpret_cast<Array3Di *>(array_); + + break; + case SHORT: + delete reinterpret_cast<Array3Ds *>(array_); + + break; + default: + assert(0); + } + } + + Image& operator= (const Image& f) { + if (this != &f) { + array_type_ = f.array_type_; + switch (array_type_) { + case BYTE: + delete reinterpret_cast<Array3Du *>(array_); + array_ = new Array3Du(*(Array3Du *)f.array_); + break; + case FLOAT: + delete reinterpret_cast<Array3Df *>(array_); + array_ = new Array3Df(*(Array3Df *)f.array_); + break; + case INT: + delete reinterpret_cast<Array3Di *>(array_); + array_ = new Array3Di(*(Array3Di *)f.array_); + break; + case SHORT: + delete reinterpret_cast<Array3Ds *>(array_); + array_ = new Array3Ds(*(Array3Ds *)f.array_); + break; + default: + assert(0); + } + } + return *this; + } + + Array3Du *AsArray3Du() const { + if (array_type_ == BYTE) { + return reinterpret_cast<Array3Du *>(array_); + } + return NULL; + } + + Array3Df *AsArray3Df() const { + if (array_type_ == FLOAT) { + return reinterpret_cast<Array3Df *>(array_); + } + return NULL; + } + + private: + DataType array_type_; + BaseArray *array_; +}; + +} // namespace libmv + +#endif // LIBMV_IMAGE_IMAGE_IMAGE_H diff --git a/intern/libmv/libmv/image/image_converter.h b/intern/libmv/libmv/image/image_converter.h new file mode 100644 index 00000000000..b3a3fa2bf8c --- /dev/null +++ b/intern/libmv/libmv/image/image_converter.h @@ -0,0 +1,81 @@ +// Copyright (c) 2009 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_IMAGE_CONVERTER_H +#define LIBMV_IMAGE_IMAGE_CONVERTER_H + +#include "libmv/image/array_nd.h" + +namespace libmv { + +// The factor comes from http://www.easyrgb.com/ +// RGB to XYZ : Y is the luminance channel +// var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722 +template<typename T> +inline T RGB2GRAY(const T r, const T g, const T b) { + return static_cast<T>(r * 0.2126 + g * 0.7152 + b * 0.0722); +} + +/* +// Specialization for the uchar type. (that do not want to work...) +template<> +inline unsigned char RGB2GRAY<unsigned char>(const unsigned char r, + const unsigned char g, + const unsigned char b) { + return (unsigned char)(r * 0.2126 + g * 0.7152 + b * 0.0722 +0.5); +}*/ + +template<class ImageIn, class ImageOut> +void Rgb2Gray(const ImageIn &imaIn, ImageOut *imaOut) { + // It is all fine to cnvert RGBA image here as well, + // all the additional channels will be nicely ignored. + assert(imaIn.Depth() >= 3); + + imaOut->Resize(imaIn.Height(), imaIn.Width(), 1); + // Convert each RGB pixel into Gray value (luminance) + + for (int j = 0; j < imaIn.Height(); ++j) { + for (int i = 0; i < imaIn.Width(); ++i) { + (*imaOut)(j, i) = RGB2GRAY(imaIn(j, i, 0) , imaIn(j, i, 1), imaIn(j, i, 2)); + } + } +} + +// Convert given float image to an unsigned char array of pixels. +template<class Image> +unsigned char *FloatImageToUCharArray(const Image &image) { + unsigned char *buffer = + new unsigned char[image.Width() * image.Height() * image.Depth()]; + + for (int y = 0; y < image.Height(); y++) { + for (int x = 0; x < image.Width(); x++) { + for (int d = 0; d < image.Depth(); d++) { + int index = (y * image.Width() + x) * image.Depth() + d; + buffer[index] = 255.0 * image(y, x, d); + } + } + } + + return buffer; +} + +} // namespace libmv + +#endif // LIBMV_IMAGE_IMAGE_CONVERTER_H diff --git a/intern/libmv/libmv/image/image_drawing.h b/intern/libmv/libmv/image/image_drawing.h new file mode 100644 index 00000000000..f50e48b75a3 --- /dev/null +++ b/intern/libmv/libmv/image/image_drawing.h @@ -0,0 +1,285 @@ +// Copyright (c) 2009 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +// Generic Image Processing Algorithm (GIPA) +// Use an ImageModel class that must implement the following : +// +// ::Contains(int y, int x) <= Tell if a point is inside or not the image +// ::operator(int y,int x) <= Modification accessor over the pixel (y,x) +// ::Width() +// ::Height() + +#ifndef LIBMV_IMAGE_IMAGE_DRAWING_H +#define LIBMV_IMAGE_IMAGE_DRAWING_H + +namespace libmv { + +/// Put the pixel in the image to the given color only if the point (xc,yc) +/// is inside the image. +template <class Image, class Color> +inline void safePutPixel(int yc, int xc, const Color & col, Image *pim) { + if (!pim) + return; + if (pim->Contains(yc, xc)) { + (*pim)(yc, xc) = col; + } +} +/// Put the pixel in the image to the given color only if the point (xc,yc) +/// is inside the image. This function support multi-channel color +/// \note The color pointer col must have size as the image depth +template <class Image, class Color> +inline void safePutPixel(int yc, int xc, const Color *col, Image *pim) { + if (!pim) + return; + if (pim->Contains(yc, xc)) { + for (int i = 0; i < pim->Depth(); ++i) + (*pim)(yc, xc, i) = *(col + i); + } +} + +// Bresenham approach to draw ellipse. +// http://raphaello.univ-fcomte.fr/ig/algorithme/ExemplesGLUt/BresenhamEllipse.htm +// Add the rotation of the ellipse. +// As the algo. use symmetry we must use 4 rotations. +template <class Image, class Color> +void DrawEllipse(int xc, int yc, int radiusA, int radiusB, + const Color &col, Image *pim, double angle = 0.0) { + int a = radiusA; + int b = radiusB; + + // Counter Clockwise rotation matrix. + double matXY[4] = { cos(angle), sin(angle), + -sin(angle), cos(angle)}; + int x, y; + double d1, d2; + x = 0; + y = b; + d1 = b*b - a*a*b + a*a/4; + + float rotX = (matXY[0] * x + matXY[1] * y); + float rotY = (matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (matXY[0] * x - matXY[1] * y); + rotY = (matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x - matXY[1] * y); + rotY = (-matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x + matXY[1] * y); + rotY = (-matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + + while (a*a*(y-.5) > b*b*(x+1)) { + if (d1 < 0) { + d1 += b*b*(2*x+3); + ++x; + } else { + d1 += b*b*(2*x+3) + a*a*(-2*y+2); + ++x; + --y; + } + rotX = (matXY[0] * x + matXY[1] * y); + rotY = (matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (matXY[0] * x - matXY[1] * y); + rotY = (matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x - matXY[1] * y); + rotY = (-matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x + matXY[1] * y); + rotY = (-matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + } + d2 = b*b*(x+.5)*(x+.5) + a*a*(y-1)*(y-1) - a*a*b*b; + while (y > 0) { + if (d2 < 0) { + d2 += b*b*(2*x+2) + a*a*(-2*y+3); + --y; + ++x; + } else { + d2 += a*a*(-2*y+3); + --y; + } + rotX = (matXY[0] * x + matXY[1] * y); + rotY = (matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (matXY[0] * x - matXY[1] * y); + rotY = (matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x - matXY[1] * y); + rotY = (-matXY[2] * x - matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + rotX = (-matXY[0] * x + matXY[1] * y); + rotY = (-matXY[2] * x + matXY[3] * y); + safePutPixel(yc + rotY, xc + rotX, col, pim); + } +} + +// Bresenham approach do not allow to draw concentric circle without holes. +// So it's better the use the Andres method. +// http://fr.wikipedia.org/wiki/Algorithme_de_tracé_de_cercle_d'Andres. +template <class Image, class Color> +void DrawCircle(int x, int y, int radius, const Color &col, Image *pim) { + Image &im = *pim; + if ( im.Contains(y + radius, x + radius) + || im.Contains(y + radius, x - radius) + || im.Contains(y - radius, x + radius) + || im.Contains(y - radius, x - radius)) { + int x1 = 0; + int y1 = radius; + int d = radius - 1; + while (y1 >= x1) { + // Draw the point for each octant. + safePutPixel( y1 + y, x1 + x, col, pim); + safePutPixel( x1 + y, y1 + x, col, pim); + safePutPixel( y1 + y, -x1 + x, col, pim); + safePutPixel( x1 + y, -y1 + x, col, pim); + safePutPixel(-y1 + y, x1 + x, col, pim); + safePutPixel(-x1 + y, y1 + x, col, pim); + safePutPixel(-y1 + y, -x1 + x, col, pim); + safePutPixel(-x1 + y, -y1 + x, col, pim); + if (d >= 2 * x1) { + d = d - 2 * x1 - 1; + x1 += 1; + } else { + if (d <= 2 * (radius - y1)) { + d = d + 2 * y1 - 1; + y1 -= 1; + } else { + d = d + 2 * (y1 - x1 - 1); + y1 -= 1; + x1 += 1; + } + } + } + } +} + +// Bresenham algorithm +template <class Image, class Color> +void DrawLine(int xa, int ya, int xb, int yb, const Color &col, Image *pim) { + Image &im = *pim; + + // If one point is outside the image + // Replace the outside point by the intersection of the line and + // the limit (either x=width or y=height). + if (!im.Contains(ya, xa) || !im.Contains(yb, xb)) { + int width = pim->Width(); + int height = pim->Height(); + const bool xdir = xa < xb, ydir = ya < yb; + float nx0 = xa, nx1 = xb, ny0 = ya, ny1 = yb, + &xleft = xdir?nx0:nx1, &yleft = xdir?ny0:ny1, + &xright = xdir?nx1:nx0, &yright = xdir?ny1:ny0, + &xup = ydir?nx0:nx1, &yup = ydir?ny0:ny1, + &xdown = ydir?nx1:nx0, &ydown = ydir?ny1:ny0; + + if (xright < 0 || xleft >= width) return; + if (xleft < 0) { + yleft -= xleft*(yright - yleft)/(xright - xleft); + xleft = 0; + } + if (xright >= width) { + yright -= (xright - width)*(yright - yleft)/(xright - xleft); + xright = width - 1; + } + if (ydown < 0 || yup >= height) return; + if (yup < 0) { + xup -= yup*(xdown - xup)/(ydown - yup); + yup = 0; + } + if (ydown >= height) { + xdown -= (ydown - height)*(xdown - xup)/(ydown - yup); + ydown = height - 1; + } + + xa = (int) xleft; + xb = (int) xright; + ya = (int) yleft; + yb = (int) yright; + } + + int xbas, xhaut, ybas, yhaut; + // Check the condition ybas < yhaut. + if (ya <= yb) { + xbas = xa; + ybas = ya; + xhaut = xb; + yhaut = yb; + } else { + xbas = xb; + ybas = yb; + xhaut = xa; + yhaut = ya; + } + // Initialize slope. + int x, y, dx, dy, incrmX, incrmY, dp, N, S; + dx = xhaut - xbas; + dy = yhaut - ybas; + if (dx > 0) { // If xhaut > xbas we will increment X. + incrmX = 1; + } else { + incrmX = -1; // else we will decrement X. + dx *= -1; + } + if (dy > 0) { // Positive slope will increment X. + incrmY = 1; + } else { // Negative slope. + incrmY = -1; + } + if (dx >= dy) { + dp = 2 * dy - dx; + S = 2 * dy; + N = 2 * (dy - dx); + y = ybas; + x = xbas; + while (x != xhaut) { + safePutPixel(y, x, col, pim); + x += incrmX; + if (dp <= 0) { // Go in direction of the South Pixel. + dp += S; + } else { // Go to the North. + dp += N; + y+=incrmY; + } + } + } else { + dp = 2 * dx - dy; + S = 2 * dx; + N = 2 * (dx - dy); + x = xbas; + y = ybas; + while (y < yhaut) { + safePutPixel(y, x, col, pim); + y += incrmY; + if (dp <= 0) { // Go in direction of the South Pixel. + dp += S; + } else { // Go to the North. + dp += N; + x += incrmX; + } + } + } + safePutPixel(y, x, col, pim); +} + +} // namespace libmv + +#endif // LIBMV_IMAGE_IMAGE_DRAWING_H diff --git a/intern/libmv/libmv/image/image_test.cc b/intern/libmv/libmv/image/image_test.cc new file mode 100644 index 00000000000..241f49f2244 --- /dev/null +++ b/intern/libmv/libmv/image/image_test.cc @@ -0,0 +1,45 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include <iostream> + +#include "libmv/image/image.h" +#include "testing/testing.h" + +using libmv::Image; +using libmv::Array3Df; + +namespace { + +TEST(Image, SimpleImageAccessors) { + Array3Df *array = new Array3Df(2, 3); + Image image(array); + EXPECT_EQ(array, image.AsArray3Df()); + EXPECT_TRUE(NULL == image.AsArray3Du()); +} + +TEST(Image, MemorySizeInBytes) { + Array3Df *array = new Array3Df(2, 3); + Image image(array); + int size = sizeof(image) + array->MemorySizeInBytes(); + EXPECT_EQ(size, image.MemorySizeInBytes()); +} + +} // namespace diff --git a/intern/libmv/libmv/image/sample.h b/intern/libmv/libmv/image/sample.h new file mode 100644 index 00000000000..24eb9ccd57d --- /dev/null +++ b/intern/libmv/libmv/image/sample.h @@ -0,0 +1,138 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_SAMPLE_H_ +#define LIBMV_IMAGE_SAMPLE_H_ + +#include "libmv/image/image.h" + +namespace libmv { + +/// Nearest neighbor interpolation. +template<typename T> +inline T SampleNearest(const Array3D<T> &image, + float y, float x, int v = 0) { + const int i = int(round(y)); + const int j = int(round(x)); + return image(i, j, v); +} + +inline void LinearInitAxis(float x, int size, + int *x1, int *x2, + float *dx) { + const int ix = static_cast<int>(x); + if (ix < 0) { + *x1 = 0; + *x2 = 0; + *dx = 1.0; + } else if (ix > size - 2) { + *x1 = size - 1; + *x2 = size - 1; + *dx = 1.0; + } else { + *x1 = ix; + *x2 = ix + 1; + *dx = *x2 - x; + } +} + +/// Linear interpolation. +template<typename T> +inline T SampleLinear(const Array3D<T> &image, float y, float x, int v = 0) { + int x1, y1, x2, y2; + float dx, dy; + + LinearInitAxis(y, image.Height(), &y1, &y2, &dy); + LinearInitAxis(x, image.Width(), &x1, &x2, &dx); + + const T im11 = image(y1, x1, v); + const T im12 = image(y1, x2, v); + const T im21 = image(y2, x1, v); + const T im22 = image(y2, x2, v); + + return T( dy * (dx * im11 + (1.0 - dx) * im12) + + (1 - dy) * (dx * im21 + (1.0 - dx) * im22)); +} + +/// Linear interpolation, of all channels. The sample is assumed to have the +/// same size as the number of channels in image. +template<typename T> +inline void SampleLinear(const Array3D<T> &image, float y, float x, T *sample) { + int x1, y1, x2, y2; + float dx, dy; + + LinearInitAxis(y, image.Height(), &y1, &y2, &dy); + LinearInitAxis(x, image.Width(), &x1, &x2, &dx); + + for (int i = 0; i < image.Depth(); ++i) { + const T im11 = image(y1, x1, i); + const T im12 = image(y1, x2, i); + const T im21 = image(y2, x1, i); + const T im22 = image(y2, x2, i); + + sample[i] = T( dy * (dx * im11 + (1.0 - dx) * im12) + + (1 - dy) * (dx * im21 + (1.0 - dx) * im22)); + } +} + +// Downsample all channels by 2. If the image has odd width or height, the last +// row or column is ignored. +// FIXME(MatthiasF): this implementation shouldn't be in an interface file +inline void DownsampleChannelsBy2(const Array3Df &in, Array3Df *out) { + int height = in.Height() / 2; + int width = in.Width() / 2; + int depth = in.Depth(); + + out->Resize(height, width, depth); + + // 2x2 box filter downsampling. + for (int r = 0; r < height; ++r) { + for (int c = 0; c < width; ++c) { + for (int k = 0; k < depth; ++k) { + (*out)(r, c, k) = (in(2 * r, 2 * c, k) + + in(2 * r + 1, 2 * c, k) + + in(2 * r, 2 * c + 1, k) + + in(2 * r + 1, 2 * c + 1, k)) / 4.0f; + } + } + } +} + +// Sample a region centered at x,y in image with size extending by half_width +// from x,y. Channels specifies the number of channels to sample from. +inline void SamplePattern(const FloatImage &image, + double x, double y, + int half_width, + int channels, + FloatImage *sampled) { + sampled->Resize(2 * half_width + 1, 2 * half_width + 1, channels); + for (int r = -half_width; r <= half_width; ++r) { + for (int c = -half_width; c <= half_width; ++c) { + for (int i = 0; i < channels; ++i) { + (*sampled)(r + half_width, c + half_width, i) = + SampleLinear(image, y + r, x + c, i); + } + } + } +} + +} // namespace libmv + +#endif // LIBMV_IMAGE_SAMPLE_H_ diff --git a/intern/libmv/libmv/image/sample_test.cc b/intern/libmv/libmv/image/sample_test.cc new file mode 100644 index 00000000000..c8a0ce470c2 --- /dev/null +++ b/intern/libmv/libmv/image/sample_test.cc @@ -0,0 +1,89 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include "libmv/image/sample.h" +#include "testing/testing.h" + +using namespace libmv; + +namespace { + +TEST(Image, Nearest) { + Array3Du image(2, 2); + image(0, 0) = 0; + image(0, 1) = 1; + image(1, 0) = 2; + image(1, 1) = 3; + EXPECT_EQ(0, SampleNearest(image, -0.4f, -0.4f)); + EXPECT_EQ(0, SampleNearest(image, 0.4f, 0.4f)); + EXPECT_EQ(3, SampleNearest(image, 0.6f, 0.6f)); + EXPECT_EQ(3, SampleNearest(image, 1.4f, 1.4f)); +} + +TEST(Image, Linear) { + Array3Df image(2, 2); + image(0, 0) = 0; + image(0, 1) = 1; + image(1, 0) = 2; + image(1, 1) = 3; + EXPECT_EQ(1.5, SampleLinear(image, 0.5, 0.5)); +} + +TEST(Image, DownsampleBy2) { + Array3Df image(2, 2); + image(0, 0) = 0; + image(0, 1) = 1; + image(1, 0) = 2; + image(1, 1) = 3; + Array3Df resampled_image; + DownsampleChannelsBy2(image, &resampled_image); + ASSERT_EQ(1, resampled_image.Height()); + ASSERT_EQ(1, resampled_image.Width()); + ASSERT_EQ(1, resampled_image.Depth()); + EXPECT_FLOAT_EQ(6./4., resampled_image(0, 0)); +} + +TEST(Image, DownsampleBy2MultiChannel) { + Array3Df image(2, 2, 3); + image(0, 0, 0) = 0; + image(0, 1, 0) = 1; + image(1, 0, 0) = 2; + image(1, 1, 0) = 3; + + image(0, 0, 1) = 5; + image(0, 1, 1) = 6; + image(1, 0, 1) = 7; + image(1, 1, 1) = 8; + + image(0, 0, 2) = 9; + image(0, 1, 2) = 10; + image(1, 0, 2) = 11; + image(1, 1, 2) = 12; + + Array3Df resampled_image; + DownsampleChannelsBy2(image, &resampled_image); + ASSERT_EQ(1, resampled_image.Height()); + ASSERT_EQ(1, resampled_image.Width()); + ASSERT_EQ(3, resampled_image.Depth()); + EXPECT_FLOAT_EQ((0+1+2+3)/4., resampled_image(0, 0, 0)); + EXPECT_FLOAT_EQ((5+6+7+8)/4., resampled_image(0, 0, 1)); + EXPECT_FLOAT_EQ((9+10+11+12)/4., resampled_image(0, 0, 2)); +} +} // namespace diff --git a/intern/libmv/libmv/image/tuple.h b/intern/libmv/libmv/image/tuple.h new file mode 100644 index 00000000000..c8dc36f2e18 --- /dev/null +++ b/intern/libmv/libmv/image/tuple.h @@ -0,0 +1,90 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#ifndef LIBMV_IMAGE_TUPLE_H +#define LIBMV_IMAGE_TUPLE_H + +#include <algorithm> + +namespace libmv { + +// A vector of elements with fixed lenght and deep copy semantics. +template <typename T, int N> +class Tuple { + public: + enum { SIZE = N }; + Tuple() {} + Tuple(T initial_value) { Reset(initial_value); } + + template <typename D> + Tuple(D *values) { Reset(values); } + + template <typename D> + Tuple(const Tuple<D, N> &b) { Reset(b); } + + template <typename D> + Tuple& operator=(const Tuple<D, N>& b) { + Reset(b); + return *this; + } + + template <typename D> + void Reset(const Tuple<D, N>& b) { Reset(b.Data()); } + + template <typename D> + void Reset(D *values) { + for (int i = 0;i < N; i++) { + data_[i] = T(values[i]); + } + } + + // Set all tuple values to the same thing. + void Reset(T value) { + for (int i = 0;i < N; i++) { + data_[i] = value; + } + } + + // Pointer to the first element. + T *Data() { return &data_[0]; } + const T *Data() const { return &data_[0]; } + + T &operator()(int i) { return data_[i]; } + const T &operator()(int i) const { return data_[i]; } + + bool operator==(const Tuple<T, N> &other) const { + for (int i = 0; i < N; ++i) { + if ((*this)(i) != other(i)) { + return false; + } + } + return true; + } + bool operator!=(const Tuple<T, N> &other) const { + return !(*this == other); + } + + private: + T data_[N]; +}; + +} // namespace libmv + +#endif // LIBMV_IMAGE_TUPLE_H diff --git a/intern/libmv/libmv/image/tuple_test.cc b/intern/libmv/libmv/image/tuple_test.cc new file mode 100644 index 00000000000..df44e5638b5 --- /dev/null +++ b/intern/libmv/libmv/image/tuple_test.cc @@ -0,0 +1,83 @@ +// Copyright (c) 2007, 2008 libmv authors. +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to +// deal in the Software without restriction, including without limitation the +// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +// sell copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +// IN THE SOFTWARE. + +#include "libmv/image/tuple.h" +#include "testing/testing.h" + +#include <algorithm> + +using libmv::Tuple; + +namespace { + +TEST(Tuple, InitConstantValue) { + Tuple<int, 3> t(5); + EXPECT_EQ(t(0), 5); + EXPECT_EQ(t(0), 5); + EXPECT_EQ(t(0), 5); +} + +TEST(Tuple, InitFromPointer) { + float vals[3] = {1.0f, 2.0f, 3.0f}; + + Tuple<float, 3> t(vals); + for (int i = 0; i < 3; i++) + EXPECT_EQ(t(i), vals[i]); + + Tuple<int, 3> b(t); + EXPECT_EQ(b(0), int(vals[0])); + EXPECT_EQ(b(1), int(vals[1])); + EXPECT_EQ(b(2), int(vals[2])); +} + +TEST(Tuple, Swap) { + unsigned char vala[3] = {1, 2, 3}; + unsigned char valb[3] = {4, 5, 6}; + + Tuple<unsigned char, 3> a(vala); + Tuple<unsigned char, 3> b(valb); + + std::swap(a, b); + + EXPECT_EQ(a(0), int(valb[0])); + EXPECT_EQ(a(1), int(valb[1])); + EXPECT_EQ(a(2), int(valb[2])); + EXPECT_EQ(b(0), int(vala[0])); + EXPECT_EQ(b(1), int(vala[1])); + EXPECT_EQ(b(2), int(vala[2])); +} + +TEST(Tuple, IsEqualOperator) { + Tuple<int, 3> a; + a(0) = 1; + a(1) = 2; + a(2) = 3; + Tuple<int, 3> b; + b(0) = 1; + b(1) = 2; + b(2) = 3; + EXPECT_TRUE(a == b); + EXPECT_FALSE(a != b); + b(1) = 5; + EXPECT_TRUE(a != b); + EXPECT_FALSE(a == b); +} + +} // namespace |