Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

Diff of /SimpleMatrix.c [000000] .. [26cac0] Maximize Restore

  Switch to unified view

a b/SimpleMatrix.c
1
2
/* autopano-sift, Automatic panorama image creation
3
 * Copyright (C) 2004 -- Sebastian Nowozin
4
 *
5
 * This program is free software released under the GNU General Public
6
 * License, which is included in this software package (doc/LICENSE).
7
 */
8
9
/* SimpleMatrix.cs
10
 *
11
 * Minimal two-dimensional matrix class implementation.
12
 *
13
 * (C) Copyright 2004 -- Sebastian Nowozin (nowozin@cs.tu-berlin.de)
14
 */
15
16
#include "AutoPanoSift.h"
17
18
SimpleMatrix* SimpleMatrix_new0 ()
19
{
20
  SimpleMatrix* self = (SimpleMatrix*)malloc(sizeof(SimpleMatrix));
21
  self->values = NULL;
22
  return self;
23
}
24
25
void SimpleMatrix_delete(SimpleMatrix* self)
26
{
27
  if (self) {
28
      DoubleMap_delete(self->values);
29
      self->values = NULL;
30
      free(self);
31
  }
32
}
33
34
void SimpleMatrix_init (SimpleMatrix* self, int yDim, int xDim)
35
{
36
  self->xDim = xDim;
37
  self->yDim = yDim;
38
  self->values = DoubleMap_new(yDim, xDim);
39
}
40
41
SimpleMatrix* SimpleMatrix_new (int yDim, int xDim)
42
{
43
  SimpleMatrix* self = SimpleMatrix_new0();
44
  self->xDim = xDim;
45
  self->yDim = yDim;
46
  self->values = DoubleMap_new(yDim, xDim);
47
  return self;
48
}
49
50
SimpleMatrix* SimpleMatrix_clone (SimpleMatrix* self)
51
{
52
  SimpleMatrix* cp = SimpleMatrix_new (self->yDim, self->xDim);
53
  
54
  int x, y;
55
  for ( y = 0 ; y < self->yDim ; ++y) {
56
      for ( x = 0 ; x < self->xDim ; ++x) {
57
      cp->values[y][x] = self->values[y][x];
58
      }
59
  }
60
  
61
  return (cp);
62
}
63
64
double SimpleMatrix_GetValue(SimpleMatrix* self, int y, int x) 
65
{
66
  return (self->values[y][x]);
67
}
68
void SimpleMatrix_SetValue(SimpleMatrix* self, int y, int x, double value )
69
{
70
  self->values[y][x] = value;
71
}
72
73
SimpleMatrix* SimpleMatrix_Mul(SimpleMatrix* m1, SimpleMatrix* m2)
74
{
75
  if (m1->xDim != m2->yDim)
76
      FatalError ("Matrixes cannot be multiplied, dimension mismatch");
77
78
  // vanilla!
79
  SimpleMatrix* res = SimpleMatrix_new (m1->yDim, m2->xDim);
80
  int y;
81
  for ( y = 0 ; y < m1->yDim ; ++y) {
82
      int x;
83
      for ( x = 0 ; x < m2->xDim ; ++x) {
84
          int k;
85
          for ( k = 0 ; k < m2->yDim ; ++k)
86
              res->values[y][x] += m1->values[y][k] * m2->values[k][x];
87
      }
88
  }
89
  
90
  return (res);
91
}
92
93
double SimpleMatrix_Dot (SimpleMatrix* self, SimpleMatrix* m)
94
{
95
  if (self->yDim != m->yDim || self->xDim != 1 || m->xDim != 1)
96
      FatalError
97
          ("Dotproduct only possible for two equal n x 1 matrices");
98
99
  double sum = 0.0;
100
101
  int y;
102
  for ( y = 0 ; y < self->yDim ; ++y)
103
      sum += self->values[y][0] * m->values[y][0];
104
105
  return (sum);
106
}
107
108
void SimpleMatrix_Negate (SimpleMatrix* self)
109
{
110
  int y;
111
  for (y = 0 ; y < self->yDim ; ++y) {
112
      int x;
113
      for ( x = 0 ; x < self->xDim ; ++x) {
114
          self->values[y][x] = -self->values[y][x];
115
      }
116
  }
117
}
118
119
void SimpleMatrix_Inverse (SimpleMatrix* self)
120
{
121
  if (self->xDim != self->yDim)
122
      FatalError("Matrix x dimension != y dimension");
123
124
  // Shipley-Coleman inversion, from
125
  // http://www.geocities.com/SiliconValley/Lab/4223/fault/ach03.html
126
  int dim = self->xDim;
127
  int i,j,k;
128
  for ( k = 0 ; k < dim ; ++k) {
129
      self->values[k][k] = - 1.0 / self->values[k][k];
130
131
      for ( i = 0 ; i < dim ; ++i) {
132
          if (i != k)
133
              self->values[i][k] *= self->values[k][k];
134
      }
135
136
      for ( i = 0 ; i < dim ; ++i) {
137
          if (i != k) {
138
              for ( j = 0 ; j < dim ; ++j) {
139
                  if (j != k)
140
                      self->values[i][j] += self->values[i][k] * self->values[k][j];
141
              }
142
          }
143
      }
144
      
145
      for ( i = 0 ; i < dim ; ++i) {
146
          if (i != k)
147
              self->values[k][i] *= self->values[k][k];
148
      }
149
      
150
  }
151
  
152
  for ( i = 0 ; i < dim ; ++i) {
153
      for ( j = 0 ; j < dim ; ++j)
154
          self->values[i][j] = -self->values[i][j];
155
  }
156
}
157
158
// The vector 'vec' is used both for input/output purposes. As input, it
159
// contains the vector v, and after this method finishes it contains x,
160
// the solution in the formula
161
//    self * x = v
162
// This matrix might get row-swapped, too.
163
void SimpleMatrix_SolveLinear (SimpleMatrix* self, SimpleMatrix* vec)
164
{
165
  if (self->xDim != self->yDim || self->yDim != vec->yDim)
166
      FatalError ("Matrix not quadratic or vector dimension mismatch");
167
  
168
  // Gaussian Elimination Algorithm, as described by
169
  // "Numerical Methods - A Software Approach", R.L. Johnston
170
  
171
  // Forward elimination with partial pivoting
172
  int x, y;
173
  for (y = 0 ; y < (self->yDim - 1) ; ++y) {
174
      
175
      // Searching for the largest pivot (to get "multipliers < 1.0 to
176
      // minimize round-off errors")
177
      int yMaxIndex = y;
178
      double yMaxValue = abs (self->values[y][y]);
179
      
180
      int py;
181
      for (py = y ; py < self->yDim ; ++py) {
182
          if (abs (self->values[py][y]) > yMaxValue) {
183
              yMaxValue = abs (self->values[py][y]);
184
              yMaxIndex = py;
185
          }
186
      }
187
188
      // if a larger row has been found, swap with the current one
189
      SimpleMatrix_SwapRow (self, y, yMaxIndex);
190
      SimpleMatrix_SwapRow (vec, y, yMaxIndex);
191
      
192
      // Now do the elimination left of the diagonal
193
      for ( py = y + 1 ; py < self->yDim ; ++py) {
194
          // always <= 1.0
195
          double elimMul = self->values[py][y] / self->values[y][y];
196
          
197
          for ( x = 0 ; x < self->xDim ; ++x)
198
              self->values[py][x] -= elimMul * self->values[y][x];
199
          
200
          // FIXME: do we really need this?
201
          vec->values[py][0] -= elimMul * vec->values[y][0];
202
      }
203
  }
204
  
205
  // Back substitution
206
  for ( y = self->yDim - 1 ; y >= 0 ; --y) {
207
      double solY = vec->values[y][0];
208
      
209
      for ( x = self->xDim - 1 ; x > y ; --x)
210
          solY -= self->values[y][x] * vec->values[x][0];
211
      
212
      vec->values[y][0] = solY / self->values[y][y];
213
  }
214
}
215
216
// Swap two rows r1, r2
217
void SimpleMatrix_SwapRow (SimpleMatrix* self, int r1, int r2)
218
{
219
  if (r1 == r2)
220
      return;
221
  
222
  int x;
223
  for ( x = 0 ; x < self->xDim ; ++x) {
224
      double temp = self->values[r1][x];
225
      self->values[r1][x] = self->values[r2][x];
226
      self->values[r2][x] = temp;
227
  }
228
}
229
230
char* SimpleMatrix_ToString (SimpleMatrix* self)
231
{
232
  int len = 5 + self->yDim * (3 + self->xDim*15);
233
  char* str = (char*)malloc(len);
234
  char* p = str;
235
236
  p += sprintf(p, "( ");
237
  int x,y;
238
  for ( y = 0 ; y < self->yDim ; ++y) {
239
      if (y > 0)
240
          p += sprintf(p, "\n  ");
241
      
242
      for ( x = 0 ; x < self->xDim ; ++x) {
243
          if (x > 0)
244
              p += sprintf(p, "  ");
245
          
246
          p += sprintf(p, "%3.015g", self->values[y][x]);
247
      }
248
  }
249
  p+= sprintf(p, " )");
250
  if (p > str + len) {
251
      FatalError("SimpleMatrix_toString overflow: len=%d, strlen=%d\n", len, p-str);
252
  }
253
254
  return str;
255
}
256
257
258
#ifdef TEST_MAIN
259
int main(int argc, char* argv[]) {
260
    SimpleMatrix * A = SimpleMatrix_new(3,3);
261
    A->values[0][0] = 1;
262
    A->values[0][1] = 2;
263
    A->values[0][2] = 3;
264
    A->values[1][0] = 4;
265
    A->values[1][1] = 5;
266
    A->values[1][2] = 6;
267
    A->values[2][0] = 7;
268
    A->values[2][1] = 9;
269
    A->values[2][2] = 8;
270
271
    WriteLine("A=%s", SimpleMatrix_ToString(A));
272
273
    SimpleMatrix * B = SimpleMatrix_new(3,1);
274
    B->values[0][0] = 1;
275
    B->values[1][0] = 2;
276
    B->values[2][0] = 3;
277
278
    WriteLine("B=%s", SimpleMatrix_ToString(B));
279
    SimpleMatrix * C = SimpleMatrix_Mul(A,B);
280
281
    WriteLine("A*B=%s", SimpleMatrix_ToString(C));
282
283
    SimpleMatrix_delete(A);
284
    SimpleMatrix_delete(B);
285
    SimpleMatrix_delete(C);
286
    return 0;
287
}
288
#endif