Rotate an image in C++ without using OpenCV functions

Question

Description :- I am trying to rotate an image without using OpenCV functions in C++. The rotation center need not be the center of the image. It could be a different point (offset from the image center). So far I followed a variety of sources to do image interpolation and I am aware of a source which does the job perfectly in MATLAB. I tried to mimic the same in C++ without OpenCV functions. But I am not getting the expected rotated image. Instead my output appears like a small horizontal line on the screen.

void RotateNearestNeighbor(cv::Mat src, double angle) {
int oldHeight = src.rows;
int oldWidth = src.cols;
int newHeight = std::sqrt(2) * oldHeight;
int newWidth = std::sqrt(2) * oldWidth;
cv::Mat output = cv::Mat(newHeight, newWidth, src.type());
double ctheta = cos(angle);
double stheta = sin(angle);

for (size_t i = 0; i < newHeight; i++) {
    for (size_t j = 0; j < newWidth; j++) {

        int oldRow = static_cast<int> ((i - newHeight / 2) * ctheta +
                                       (j - newWidth / 2) * stheta + oldHeight / 2);
        int oldCol = static_cast<int> (-(i - newHeight / 2) * stheta +
                                       (j - newWidth / 2) * ctheta + oldWidth / 2);

        if (oldRow > 0 && oldCol > 0 && oldRow <= oldHeight && oldCol <= oldWidth)
            output.at<cv::Vec3b>(i, j) = src.at<cv::Vec3b>(oldRow, oldCol);
        else
            output.at<cv::Vec3b>(i, j) = cv::Vec3b(0, 0, 0);
    }
}
cv::imshow("Rotated cat", output);
}

The following are my input (left side) and output (right side) images

---

UPDATE : -

After being inspired by many answers related to this question and also the most elaborate, helpful and generous answer below, I could fix my OpenCV code to get the desired result.

Modified Code :

// Trivial constant
constexpr double Pi = 3.1415926535897932384626433832795;

/*!
* rief Function to generate transformation matrix
* param angle is the angle of rotation from user input
* param pivot is the amount of translation in x and y axes
* 
eturn translation matrix
*/
cv::Mat CreateTransMat(double angle, std::pair<int, int> &pivot) {
    angle = Pi * angle / 180;
    return (cv::Mat_<double>(3, 3) << cos(angle), -sin(angle), pivot.first,
            sin(angle), cos(angle), pivot.second, 0, 0, 1);
}

/*!
* rief Function to apply coordinate transform from destination to     source
* param inv_mat being the inverse transformation matrix for the transform needed
* 
eturn pos being the homogeneous coordinates for transformation
*/
cv::Mat CoordTransform(const cv::Mat &inv_mat, const cv::Mat &pos) {
    assert(inv_mat.cols == pos.rows);
    cv::Mat trans_mat = inv_mat * pos;
    return (cv::Mat_<double>(1, 2) <<
            trans_mat.at<double>(0, 0) / trans_mat.at<double>(0, 2),
            trans_mat.at<double>(0, 1) / trans_mat.at<double>(0, 2));
}

/*!
* rief Function to transform an image based on a rotation angle and translation
         matrix. When rotation and translation happen at the same time, the
         two matrices can be combined
* param src being source image
* param dest being destination image
* param trans_mat being the transformation (rotation/ translation) matrix
*/
void ImageTransform(const cv::Mat &src, const cv::Mat &trans_mat, cv::Mat &dest) {
    int src_rows = src.rows;
    int src_cols = src.cols;
    int dest_rows = dest.rows;
    int dest_cols = dest.cols;
    const cv::Mat inverse_mat = trans_mat.inv();
    //#pragma omp parallel for simd
    for (int row = 0; row < dest_rows; row++) {
        //#pragma omp parallel for simd
        for (int col = 0; col < dest_cols; col++) {
            cv::Mat src_pos = CoordTransform(inverse_mat,
                                         (cv::Mat_<double>(3, 1) << col, row, 1));
            const int x_actual = static_cast<int>(src_pos.at<double>(0, 0) + 0.5);
            const int y_actual = static_cast<int>(src_pos.at<double>(0, 1) + 0.5);

            if (x_actual >= 0 && x_actual < src_cols &&
                y_actual >= 0 && y_actual < src_rows)
                dest.at<cv::Vec3b>(row, col) = src.at<cv::Vec3b>(y_actual, x_actual);
            else
                dest.at<cv::Vec3b>(row, col) = cv::Vec3b(0, 0, 0);
        }
    }    
}

/*!
* rief User manual for command-line args input
*/
void Usage() {
    std::cout << "COMMAND INPUT : - 

" <<
              "          ./ImageTransform <image> <rotation-angle>" <<
              std::endl;
}
/*!
* rief main function to read a user input location for an image and then apply the
         required transformations (rotation / translation)
*/
int main(int argc, char *argv[])
{
    auto start = std::chrono::steady_clock::now();
    if (argc == 0 || argc < 3)
        Usage();
    else {
        double degree = std::stod(argv[2]);
        double angle = degree * CV_PI / 180.;
        cv::Mat src_img = cv::imread(argv[1]);
        std::pair<int, int> null_trans = std::make_pair(0, 0);
        std::pair<int, int> translation_initial =
            std::make_pair(src_img.cols / 2 + 1, src_img.rows / 2 + 1);
        std::pair<int, int> translation_final =
            std::make_pair(0, -src_img.rows / 2 - 4);
        if (!src_img.data)
        {
            std::cout << "image null" << std::endl;
            cv::waitKey(0);
        }
        cv::imshow("Source", src_img);
        cv::Mat dest_img = cv::Mat(static_cast<int>(2 * src_img.rows),
                                   static_cast<int>(2 * src_img.cols),
                                   src_img.type());
        cv::Mat trans_mat1 = CreateTransMat(degree, translation_initial);
        ImageTransform(src_img, trans_mat1, dest_img);
        cv::imshow("Interim", dest_img);
        cv::Mat interim_img = dest_img;
        dest_img.release();
        dest_img = cv::Mat(src_img.rows, src_img.cols, src_img.type());
        cv::Mat trans_mat2 = CreateTransMat(0, translation_final);
        ImageTransform(interim_img, trans_mat2, dest_img);
        cv::imshow("Final image", dest_img);
        cv::waitKey(10);
    }
    auto end = std::chrono::steady_clock::now();
    auto diff = end - start;
    std::cout << std::chrono::duration <double, std::milli> (diff).count() <<
              " ms" << std::endl;
}

Input image

Rotated image

Answer 1

First, I have to admit I agree with generic_opto_guy:

The approach with the loop looks good, so we would need to check the math. On thing I noticed: if (oldRow > 0 && oldCol > 0 && oldRow <= oldHeight && oldCol <= oldWidth) implies you start indexing with 1. I belife that opencv starts indexing with 0.

For all that, I couldn't resist to answer. (May be, it's just an image phase of mine.)

Instead of fiddling with sin() and cos(), I would recommend to use matrix transformation. At the first glance, this might appear over-engineered but later you will recognize that it bears much more flexibility. With a transformation matrix, you can express a lot of transformations (translation, rotation, scaling, shearing, projection) as well as combining multiple transformations into one matrix.

(A teaser for what is possible: SO: How to paint / deform a QImage in 2D?)

In an image, the pixels may be addressed by 2d coordinates. Hence a 2×2 matrix comes into mind but a 2×2 matrix cannot express translations. For this, homogeneous coordinates has been introduced – a math trick to handle positions and directions in the same space by extending the dimension by one.

To make it short, a 2d position (x, y) has the homogeneous coordinates (x, y, 1).

A position transformed with a transformation matrix:

v′ = M · v.

This may or may not change the value of third component. To convert the homogeneous coordinate to 2D position again, x and y has to be divided by 3^rd component.

Vec2 transform(const Mat3x3 &mat, const Vec2 &pos)
{
  const Vec3 pos_ = mat * Vec3(pos, 1.0);
  return Vec2(pos_.x / pos_.z, pos_.y / pos_.z);
}

To transform a source image into a destination image, the following function can be used:

void transform(
  const Image &imgSrc, const Mat3x3 &mat, Image &imgDst,
  int rgbFail = 0x808080)
{
  const Mat3x3 matInv = invert(mat);
  for (int y = 0; y < imgDst.h(); ++y) {
    for (int x = 0; x < imgDst.w(); ++x) {
      const Vec2 pos = transform(matInv, Vec2(x, y));
      const int xSrc = (int)(pos.x + 0.5), ySrc = (int)(pos.y + 0.5);
      imgDst.setPixel(x, y,
        xSrc >= 0 && xSrc < imgSrc.w() && ySrc >= 0 && ySrc < imgSrc.h()
        ? imgSrc.getPixel(xSrc, ySrc)
        : rgbFail);
    }
  }
}

Note:

The transformation matrix mat describes the transformation of source image coordinates to destination image coordinates. The nested loops iterate over destination image. Hence, the inverse matrix (representing the reverse transformation) has to be used to get the corresponding source image coordinates which map to the current destination coordinates.

… and the matrix constructor for the rotation:

enum ArgInitRot { InitRot };

template <typename VALUE>
struct Mat3x3T {
  union {
    VALUE comp[3 * 3];
    struct {
      VALUE _00, _01, _02;
      VALUE _10, _11, _12;
      VALUE _20, _21, _22;
    };
  };

  // constructor to build a matrix for rotation
  Mat3x3T(ArgInitRot, VALUE angle):
    _00(std::cos(angle)), _01(-std::sin(angle)), _02((VALUE)0),
    _10(std::sin(angle)), _11( std::cos(angle)), _12((VALUE)0),
    _20(       (VALUE)0), _21(        (VALUE)0), _22((VALUE)1)
  { }

can be used to construct a rotation with angle (in degree):

Mat3x3T<double> mat(InitRot, degToRad(30.0));

Note:

I would like to emphasize how the transformed coordinates are used:

      const Vec2 pos = transform(matInv, Vec2(x, y));
      const int xSrc = (int)(pos.x + 0.5), ySrc = (int)(pos.y + 0.5);

Rounding the results to yield one discrete pixel position is actually what is called Nearest Neighbour. Alternatively, the now discarded fractional parts could be used for a linear interpolation between neighbour pixels.

---

To make a small sample, I first copied image.h, image.cc, imagePPM.h, and imagePPM.cc from another answer I wrote recently. (The PPM file format has been used as it needs minimal code for file I/O.)

Next, I used linMath.h (my minimal math collection for 3D transformations) to make a minimal math collection for 2D transformations – linMath.h:

#ifndef LIN_MATH_H
#define LIN_MATH_H

#include <iostream>
#include <cassert>
#include <cmath>

extern const double Pi;

template <typename VALUE>
inline VALUE degToRad(VALUE angle)
{
  return (VALUE)Pi * angle / (VALUE)180;
}

template <typename VALUE>
inline VALUE radToDeg(VALUE angle)
{
  return (VALUE)180 * angle / (VALUE)Pi;
}

enum ArgNull { Null };

template <typename VALUE>
struct Vec2T {
  typedef VALUE Value;
  Value x, y;
  // default constructor (leaving elements uninitialized)
  Vec2T() { }
  Vec2T(ArgNull): x((Value)0), y((Value)0) { }
  Vec2T(Value x, Value y): x(x), y(y) { }
};

typedef Vec2T<float> Vec2f;
typedef Vec2T<double> Vec2;

template <typename VALUE>
struct Vec3T {
  typedef VALUE Value;
  Value x, y, z;
  // default constructor (leaving elements uninitialized)
  Vec3T() { }
  Vec3T(ArgNull): x((Value)0), y((Value)0), z((Value)0) { }
  Vec3T(Value x, Value y, Value z): x(x), y(y), z(z) { }
  Vec3T(const Vec2T<Value> &xy, Value z): x(xy.x), y(xy.y), z(z) { }
  explicit operator Vec2T<Value>() const { return Vec2T<Value>(x, y); }
  const Vec2f xy() const { return Vec2f(x, y); }
  const Vec2f xz() const { return Vec2f(x, z); }
  const Vec2f yz() const { return Vec2f(y, z); }
};

typedef Vec3T<float> Vec3f;
typedef Vec3T<double> Vec3;

enum ArgInitIdent { InitIdent };
enum ArgInitTrans { InitTrans };
enum ArgInitRot { InitRot };
enum ArgInitScale { InitScale };
enum ArgInitFrame { InitFrame };

template <typename VALUE>
struct Mat3x3T {
  union {
    VALUE comp[3 * 3];
    struct {
      VALUE _00, _01, _02;
      VALUE _10, _11, _12;
      VALUE _20, _21, _22;
    };
  };

  // default constructor (leaving elements uninitialized)
  Mat3x3T() { }
  // constructor to build a matrix by elements
  Mat3x3T(
    VALUE _00, VALUE _01, VALUE _02,
    VALUE _10, VALUE _11, VALUE _12,
    VALUE _20, VALUE _21, VALUE _22):
    _00(_00), _01(_01), _02(_02),
    _10(_10), _11(_11), _12(_12),
    _20(_20), _21(_21), _22(_22)
  { }
  // constructor to build an identity matrix
  Mat3x3T(ArgInitIdent):
    _00((VALUE)1), _01((VALUE)0), _02((VALUE)0),
    _10((VALUE)0), _11((VALUE)1), _12((VALUE)0),
    _20((VALUE)0), _21((VALUE)0), _22((VALUE)1)
  { }
  // constructor to build a matrix for translation
  Mat3x3T(ArgInitTrans, const Vec2T<VALUE> &t):
    _00((VALUE)1), _01((VALUE)0), _02((VALUE)t.x),
    _10((VALUE)0), _11((VALUE)1), _12((VALUE)t.y),
    _20((VALUE)0), _21((VALUE)0), _22((VALUE)1)
  { }
  // constructor to build a matrix for rotation
  Mat3x3T(ArgInitRot, VALUE angle):
    _00(std::cos(angle)), _01(-std::sin(angle)), _02((VALUE)0),
    _10(std::sin(angle)), _11( std::cos(angle)), _12((VALUE)0),
    _20(       (VALUE)0), _21(        (VALUE)0), _22((VALUE)1)
  { }
  // constructor to build a matrix for translation/rotation
  Mat3x3T(ArgInitFrame, const Vec2T<VALUE> &t, VALUE angle):
    _00(std::cos(angle)), _01(-std::sin(angle)), _02((VALUE)t.x),
    _10(std::sin(angle)), _11( std::cos(angle)), _12((VALUE)t.y),
    _20(       (VALUE)0), _21(        (VALUE)0), _22((VALUE)1)
  { }
  // constructor to build a matrix for scaling
  Mat3x3T(ArgInitScale, VALUE sx, VALUE sy):
    _00((VALUE)sx), _01( (VALUE)0), _02((VALUE)0),
    _10( (VALUE)0), _11((VALUE)sy), _12((VALUE)0),
    _20( (VALUE)0), _21( (VALUE)0), _22((VALUE)1)
  { }
  // operator to allow access with [][]
  VALUE* operator [] (int i)
  {
    assert(i >= 0 && i < 3);
    return comp + 3 * i;
  }
  // operator to allow access with [][]
  const VALUE* operator [] (int i) const
  {
    assert(i >= 0 && i < 3);
    return comp + 3 * i;
  }

  // multiply matrix with matrix -> matrix
  Mat3x3T operator * (const Mat3x3T &mat) const
  {
    return Mat3x3T(
      _00 * mat._00 + _01 * mat._10 + _02 * mat._20,
      _00 * mat._01 + _01 * mat._11 + _02 * mat._21,
      _00 * mat._02 + _01 * mat._12 + _02 * mat._22,
      _10 * mat._00 + _11 * mat._10 + _12 * mat._20,
      _10 * mat._01 + _11 * mat._11 + _12 * mat._21,
      _10 * mat._02 + _11 * mat._12 + _12 * mat._22,
      _20 * mat._00 + _21 * mat._10 + _22 * mat._20,
      _20 * mat._01 + _21 * mat._11 + _22 * mat._21,
      _20 * mat._02 + _21 * mat._12 + _22 * mat._22);
  }
  // multiply matrix with vector -> vector
  Vec3T<VALUE> operator * (const Vec3T<VALUE> &vec) const
  {
    return Vec3T<VALUE>(
      _00 * vec.x + _01 * vec.y + _02 * vec.z,
      _10 * vec.x + _11 * vec.y + _12 * vec.z,
      _20 * vec.x + _21 * vec.y + _22 * vec.z);
  }
};

typedef Mat3x3T<float> Mat3x3f;
typedef Mat3x3T<double> Mat3x3;

template <typename VALUE>
std::ostream& operator<<(std::ostream &out, const Mat3x3T<VALUE> &m)
{
  return out
    << m._00 << '' << m._01 << '' << m._02 << '
'
    << m._10 << '' << m._11 << '' << m._12 << '
'
    << m._20 << '' << m._21 << '' << m._22 << '
';
}

/* computes determinant of a matrix.
 *
 * det = |M|
 *
 * mat ... the matrix
 */
template <typename VALUE>
VALUE determinant(const Mat3x3T<VALUE> &mat)
{
  return mat._00 * mat._11 * mat._22
    + mat._01 * mat._12 * mat._20
    + mat._02 * mat._10 * mat._21
    - mat._20 * mat._11 * mat._02
    - mat._21 * mat._12 * mat._00
    - mat._22 * mat._10 * mat._01;
}

/* returns the inverse of a regular matrix.
 *
 * mat matrix to invert
 * eps epsilon for regularity of matrix
 */
template <typename VALUE>
Mat3x3T<VALUE> invert(
  const Mat3x3T<VALUE> &mat, VALUE eps = (VALUE)1E-10)
{
  assert(eps >= (VALUE)0);
  // compute determinant and check that it its unequal to 0
  // (Otherwise, matrix is singular!)
  const VALUE det = determinant(mat);
  if (std::abs(det) < eps) throw std::domain_error("Singular matrix!");
  // reciproke of determinant
  const VALUE detInvPos = (VALUE)1 / det, detInvNeg = -detInvPos;
  // compute each element by determinant of sub-matrix which is build
  // striking out row and column of pivot element itself
  // BTW, the determinant is multiplied with -1 when sum of row and column
  // index is odd (chess board rule)
  // (This is usually called cofactor of related element.)
  // transpose matrix and multiply with 1/determinant of original matrix
  return Mat3x3T<VALUE>(
    detInvPos * (mat._11 * mat._22 - mat._12 * mat._21),
    detInvNeg * (mat._