matrix.cpp 4.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205
  1. /*
  2. * DoRayMe - a quick and dirty Raytracer
  3. * Matrix implementation
  4. *
  5. * Created by Manoël Trapier
  6. * Copyright (c) 2020 986-Studio.
  7. *
  8. */
  9. #include <stdio.h>
  10. #include <matrix.h>
  11. #include <tuple.h>
  12. #include <math_helper.h>
  13. Matrix::Matrix(int width)
  14. {
  15. int i;
  16. this->size = width;
  17. for(i = 0; i < width*width; i++)
  18. {
  19. this->data[i] = 0;
  20. }
  21. };
  22. Matrix::Matrix(double values[], int width)
  23. {
  24. int x, y;
  25. this->size = width;
  26. for(y = 0; y < this->size; y++)
  27. {
  28. for (x = 0 ; x < this->size ; x++)
  29. {
  30. this->data[this->size * x + y] = values[this->size * x + y];
  31. }
  32. }
  33. };
  34. bool Matrix::operator==(const Matrix &b) const
  35. {
  36. int i;
  37. if (this->size != b.size)
  38. {
  39. /* If they are not the same size don't even bother */
  40. return false;
  41. }
  42. for(i = 0; i < this->size*this->size; i++)
  43. {
  44. if (!double_equal(this->data[i], b.data[i]))
  45. {
  46. return false;
  47. }
  48. }
  49. return true;
  50. }
  51. bool Matrix::operator!=(const Matrix &b) const
  52. {
  53. int i;
  54. if (this->size != b.size)
  55. {
  56. /* If they are not the same size don't even bother */
  57. return true;
  58. }
  59. for(i = 0; i < this->size*this->size; i++)
  60. {
  61. if (!double_equal(this->data[i], b.data[i]))
  62. {
  63. return true;
  64. }
  65. }
  66. return false;
  67. }
  68. Matrix Matrix::operator*(const Matrix &b) const
  69. {
  70. int x, y, k;
  71. Matrix ret = Matrix(this->size);
  72. if (this->size == b.size)
  73. {
  74. for (y = 0 ; y < this->size ; y++)
  75. {
  76. for (x = 0 ; x < this->size ; x++)
  77. {
  78. double v = 0;
  79. for (k = 0 ; k < this->size ; k++)
  80. {
  81. v += this->get(x, k) * b.get(k, y);
  82. }
  83. ret.set(x, y, v);
  84. }
  85. }
  86. }
  87. return ret;
  88. }
  89. /* TODO: Check if we can optimise this function. It is called a lot */
  90. Tuple Matrix::operator*(const Tuple &b) const
  91. {
  92. return Tuple(b.x * this->get(0, 0) + b.y * this->get(0, 1) + b.z * this->get(0, 2) + b.w * this->get(0, 3),
  93. b.x * this->get(1, 0) + b.y * this->get(1, 1) + b.z * this->get(1, 2) + b.w * this->get(1, 3),
  94. b.x * this->get(2, 0) + b.y * this->get(2, 1) + b.z * this->get(2, 2) + b.w * this->get(2, 3),
  95. b.x * this->get(3, 0) + b.y * this->get(3, 1) + b.z * this->get(3, 2) + b.w * this->get(3, 3));
  96. }
  97. Matrix Matrix::identity()
  98. {
  99. int i;
  100. for(i = 0; i < this->size; i++)
  101. {
  102. this->set(i, i, 1);
  103. }
  104. return *this;
  105. }
  106. Matrix Matrix::transpose()
  107. {
  108. int x, y;
  109. Matrix ret = Matrix(this->size);
  110. #pragma omp parallel for simd private(y, x)
  111. for (y = 0 ; y < this->size ; y++)
  112. {
  113. for (x = 0 ; x < this->size ; x++)
  114. {
  115. ret.set(y, x, this->get(x, y));
  116. }
  117. }
  118. return ret;
  119. }
  120. Matrix Matrix::submatrix(int row, int column)
  121. {
  122. int i, j;
  123. int x = 0, y = 0;
  124. Matrix ret = Matrix(this->size - 1);
  125. for (i = 0 ; i < this->size ; i++)
  126. {
  127. if (i == row)
  128. {
  129. /* Skip that row */
  130. continue;
  131. }
  132. y = 0;
  133. for (j = 0 ; j < this->size ; j++)
  134. {
  135. if (j == column)
  136. {
  137. /* skip that column */
  138. continue;
  139. }
  140. ret.set(x, y, this->get(i, j));
  141. y++;
  142. }
  143. x++;
  144. }
  145. return ret;
  146. }
  147. double Matrix::determinant()
  148. {
  149. double det = 0;
  150. if (this->size == 2)
  151. {
  152. det = this->data[0] * this->data[3] - this->data[1] * this->data[2];
  153. }
  154. else
  155. {
  156. int col;
  157. for(col = 0; col < this->size; col++)
  158. {
  159. det = det + this->get(0, col) * this->cofactor(0, col);
  160. }
  161. }
  162. return det;
  163. }
  164. Matrix Matrix::inverse()
  165. {
  166. Matrix ret = Matrix(this->size);
  167. if (this->isInvertible())
  168. {
  169. int row, col;
  170. double c;
  171. for (row = 0; row < this->size; row++)
  172. {
  173. for (col = 0; col < this->size; col++)
  174. {
  175. c = this->cofactor(row, col);
  176. ret.set(col, row, c / this->determinant());
  177. }
  178. }
  179. }
  180. return ret;
  181. }