← Back to branch summary

~ubuntu-branches/ubuntu/vivid/rawstudio/vivid

« back to all changes in this revision

Viewing changes to librawstudio/rs-math.c

  • Committer: Bazaar Package Importer
  • Author(s): Bernd Zeimetz
  • Date: 2011-07-28 17:36:32 UTC
  • mfrom: (2.1.11 upstream)
  • Revision ID: james.westby@ubuntu.com-20110728173632-5czluz9ye3c83zc5
Tags: 2.0-1
* [3750b2cf] Merge commit 'upstream/2.0'
* [63637468] Removing Patch, not necessary anymore.
* [2fb580dc] Add new build-dependencies.
* [c57d953b] Run dh_autoreconf due to patches in configure.in
* [13febe39] Add patch to remove the libssl requirement.
* [5ae773fe] Replace libjpeg62-dev by libjpeg8-dev :)
* [1969d755] Don't build static libraries.
* [7cfe0a2e] Add a patch to fix the plugin directory path.
  As plugins are shared libraries, they need to go into /usr/lib,
  not into /usr/share.
  Thanks to Andrew McMillan
* [c1d0d9dd] Don't install .la files for all plugins and libraries.

Show diffs side-by-side

added added

removed removed

Lines of Context:
librawstudio/rs-math.c
 
1
/*
 
2
 * * Copyright (C) 2006-2011 Anders Brander <anders@brander.dk>, 
 
3
 * * Anders Kvist <akv@lnxbx.dk> and Klaus Post <klauspost@gmail.com>
 
4
 *
 
5
 * This program is free software; you can redistribute it and/or
 
6
 * modify it under the terms of the GNU General Public License
 
7
 * as published by the Free Software Foundation; either version 2
 
8
 * of the License, or (at your option) any later version.
 
9
 *
 
10
 * This program is distributed in the hope that it will be useful,
 
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 * GNU General Public License for more details.
 
14
 *
 
15
 * You should have received a copy of the GNU General Public License
 
16
 * along with this program; if not, write to the Free Software
 
17
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
18
 */
 
19
 
 
20
#include <stdio.h>
 
21
#include <math.h>
 
22
#include <stdlib.h>
 
23
#include "rs-math.h"
 
24
 
 
25
/* luminance weights, notice that these is used for linear data */
 
26
#define RLUM (0.3086)
 
27
#define GLUM (0.6094)
 
28
#define BLUM (0.0820)
 
29
 
 
30
static void matrix4_zshear (RS_MATRIX4 *matrix, double dx, double dy);
 
31
static void matrix4_xrotate(RS_MATRIX4 *matrix, double rs, double rc);
 
32
static void matrix4_yrotate(RS_MATRIX4 *matrix, double rs, double rc);
 
33
static void matrix4_zrotate(RS_MATRIX4 *matrix, double rs, double rc);
 
34
static void matrix4_affine_transform_3dpoint(RS_MATRIX4 *matrix, double x, double y, double z, double *tx, double *ty, double *tz);
 
35
 
 
36
void
 
37
printvec(const char *str, const RS_VECTOR3 *vec)
 
38
{
 
39
        printf("%s: [ %.05f %.05f %.05f ]\n", str, vec->x, vec->y, vec->z);
 
40
}
 
41
 
 
42
void
 
43
printmat3(RS_MATRIX3 *mat)
 
44
{
 
45
        int y;
 
46
 
 
47
        printf("M: matrix(\n");
 
48
        for(y=0; y<3; y++)
 
49
        {
 
50
                printf("\t[ %f, ",mat->coeff[y][0]);
 
51
                printf("%f, ",mat->coeff[y][1]);
 
52
                printf("%f ",mat->coeff[y][2]);
 
53
                printf("],\n");
 
54
        }
 
55
        printf(")\n");
 
56
}
 
57
 
 
58
void
 
59
printmat(RS_MATRIX4 *mat)
 
60
{
 
61
        int x, y;
 
62
 
 
63
        for(y=0; y<4; y++)
 
64
        {
 
65
                for(x=0; x<4; x++)
 
66
                        printf("%f ",mat->coeff[y][x]);
 
67
                printf("\n");
 
68
        }
 
69
        printf("\n");
 
70
}
 
71
 
 
72
float
 
73
vector3_max(const RS_VECTOR3 *vec)
 
74
{
 
75
        float max = 0.0;
 
76
 
 
77
        max = MAX(max, vec->x);
 
78
        max = MAX(max, vec->y);
 
79
        max = MAX(max, vec->z);
 
80
 
 
81
        return max;
 
82
}
 
83
 
 
84
RS_MATRIX3
 
85
vector3_as_diagonal(const RS_VECTOR3 *vec)
 
86
{
 
87
        RS_MATRIX3 result;
 
88
        matrix3_identity(&result);
 
89
        result.coeff[0][0] = vec->x;
 
90
        result.coeff[1][1] = vec->y;
 
91
        result.coeff[2][2] = vec->z;
 
92
 
 
93
        return result;
 
94
}
 
95
 
 
96
void
 
97
matrix4_identity (RS_MATRIX4 *matrix)
 
98
{
 
99
  static const RS_MATRIX4 identity = { { { 1.0, 0.0, 0.0, 0.0 },
 
100
                                          { 0.0, 1.0, 0.0, 0.0 },
 
101
                                          { 0.0, 0.0, 1.0, 0.0 },
 
102
                                          { 0.0, 0.0, 0.0, 1.0 } } };
 
103
 
 
104
  *matrix = identity;
 
105
}
 
106
 
 
107
void
 
108
matrix4_multiply(const RS_MATRIX4 *left, RS_MATRIX4 *right, RS_MATRIX4 *result)
 
109
{
 
110
  int i, j;
 
111
  RS_MATRIX4 tmp;
 
112
  double t1, t2, t3, t4;
 
113
 
 
114
  for (i = 0; i < 4; i++)
 
115
    {
 
116
      t1 = left->coeff[i][0];
 
117
      t2 = left->coeff[i][1];
 
118
      t3 = left->coeff[i][2];
 
119
      t4 = left->coeff[i][3];
 
120
 
 
121
      for (j = 0; j < 4; j++)
 
122
        {
 
123
          tmp.coeff[i][j]  = t1 * right->coeff[0][j];
 
124
          tmp.coeff[i][j] += t2 * right->coeff[1][j];
 
125
          tmp.coeff[i][j] += t3 * right->coeff[2][j];
 
126
          tmp.coeff[i][j] += t4 * right->coeff[3][j];
 
127
        }
 
128
    }
 
129
  *result = tmp;
 
130
}
 
131
 
 
132
/* copied almost verbatim from  dcraw.c:pseudoinverse()
 
133
 - but this one doesn't transpose */
 
134
void
 
135
matrix4_color_invert(const RS_MATRIX4 *in, RS_MATRIX4 *out)
 
136
{
 
137
        RS_MATRIX4 tmp;
 
138
        double work[3][6], num;
 
139
        int i, j, k;
 
140
 
 
141
        matrix4_identity(&tmp);
 
142
 
 
143
        for (i=0; i < 3; i++)
 
144
        {
 
145
                for (j=0; j < 6; j++)
 
146
                        work[i][j] = j == i+3;
 
147
                for (j=0; j < 3; j++)
 
148
                        for (k=0; k < 3; k++)
 
149
                                work[i][j] += in->coeff[k][i] * in->coeff[k][j];
 
150
        }
 
151
        for (i=0; i < 3; i++)
 
152
        {
 
153
                num = work[i][i];
 
154
                for (j=0; j < 6; j++)
 
155
                        work[i][j] /= num;
 
156
                for (k=0; k < 3; k++)
 
157
                {
 
158
                        if (k==i)
 
159
                                continue;
 
160
                        num = work[k][i];
 
161
                        for (j=0; j < 6; j++)
 
162
                                work[k][j] -= work[i][j] * num;
 
163
                }
 
164
        }
 
165
        for (i=0; i < 3; i++)
 
166
                for (j=0; j < 3; j++)
 
167
                        for (tmp.coeff[i][j]=k=0; k < 3; k++)
 
168
                                tmp.coeff[i][j] += work[j][k+3] * in->coeff[i][k];
 
169
        for (i=0; i < 4; i++)
 
170
                for (j=0; j < 4; j++)
 
171
                        out->coeff[i][j] = tmp.coeff[j][i];
 
172
}
 
173
 
 
174
static void
 
175
matrix4_zshear (RS_MATRIX4 *matrix, double dx, double dy)
 
176
{
 
177
  RS_MATRIX4 zshear;
 
178
 
 
179
  zshear.coeff[0][0] = 1.0;
 
180
  zshear.coeff[1][0] = 0.0;
 
181
  zshear.coeff[2][0] = dx;
 
182
  zshear.coeff[3][0] = 0.0;
 
183
 
 
184
  zshear.coeff[0][1] = 0.0;
 
185
  zshear.coeff[1][1] = 1.0;
 
186
  zshear.coeff[2][1] = dy;
 
187
  zshear.coeff[3][1] = 0.0;
 
188
 
 
189
  zshear.coeff[0][2] = 0.0;
 
190
  zshear.coeff[1][2] = 0.0;
 
191
  zshear.coeff[2][2] = 1.0;
 
192
  zshear.coeff[3][2] = 0.0;
 
193
 
 
194
  zshear.coeff[0][3] = 0.0;
 
195
  zshear.coeff[1][3] = 0.0;
 
196
  zshear.coeff[2][3] = 0.0;
 
197
  zshear.coeff[3][3] = 1.0;
 
198
 
 
199
  matrix4_multiply(&zshear, matrix, matrix);
 
200
}
 
201
 
 
202
void
 
203
matrix4_to_matrix4int(RS_MATRIX4 *matrix, RS_MATRIX4Int *matrixi)
 
204
{
 
205
        int a,b;
 
206
        for(a=0;a<4;a++)
 
207
                for(b=0;b<4;b++)
 
208
                {
 
209
                        /* Check that a (unsigned 16 bit) x (matrix coeff shifted up RESOLUTION) fits within signed 32 bit value */
 
210
                        /* Adjust this if MATRIX_RESOLUTION is adjusted */
 
211
                        g_assert((matrix->coeff[a][b] < 16.0) && (matrix->coeff[a][b] > -16.0));
 
212
 
 
213
                        matrixi->coeff[a][b] = (int) (matrix->coeff[a][b] * (double) (1<<MATRIX_RESOLUTION));
 
214
                }
 
215
}
 
216
 
 
217
static void
 
218
matrix4_xrotate(RS_MATRIX4 *matrix, double rs, double rc)
 
219
{
 
220
  RS_MATRIX4 tmp;
 
221
 
 
222
  tmp.coeff[0][0] = 1.0;
 
223
  tmp.coeff[1][0] = 0.0;
 
224
  tmp.coeff[2][0] = 0.0;
 
225
  tmp.coeff[3][0] = 0.0;
 
226
 
 
227
  tmp.coeff[0][1] = 0.0;
 
228
  tmp.coeff[1][1] = rc;
 
229
  tmp.coeff[2][1] = rs;
 
230
  tmp.coeff[3][1] = 0.0;
 
231
 
 
232
  tmp.coeff[0][2] = 0.0;
 
233
  tmp.coeff[1][2] = -rs;
 
234
  tmp.coeff[2][2] = rc;
 
235
  tmp.coeff[3][2] = 0.0;
 
236
 
 
237
  tmp.coeff[0][3] = 0.0;
 
238
  tmp.coeff[1][3] = 0.0;
 
239
  tmp.coeff[2][3] = 0.0;
 
240
  tmp.coeff[3][3] = 1.0;
 
241
 
 
242
  matrix4_multiply(&tmp, matrix, matrix);
 
243
}
 
244
 
 
245
static void
 
246
matrix4_yrotate(RS_MATRIX4 *matrix, double rs, double rc)
 
247
{
 
248
  RS_MATRIX4 tmp;
 
249
 
 
250
  tmp.coeff[0][0] = rc;
 
251
  tmp.coeff[1][0] = 0.0;
 
252
  tmp.coeff[2][0] = -rs;
 
253
  tmp.coeff[3][0] = 0.0;
 
254
 
 
255
  tmp.coeff[0][1] = 0.0;
 
256
  tmp.coeff[1][1] = 1.0;
 
257
  tmp.coeff[2][1] = 0.0;
 
258
  tmp.coeff[3][1] = 0.0;
 
259
 
 
260
  tmp.coeff[0][2] = rs;
 
261
  tmp.coeff[1][2] = 0.0;
 
262
  tmp.coeff[2][2] = rc;
 
263
  tmp.coeff[3][2] = 0.0;
 
264
 
 
265
  tmp.coeff[0][3] = 0.0;
 
266
  tmp.coeff[1][3] = 0.0;
 
267
  tmp.coeff[2][3] = 0.0;
 
268
  tmp.coeff[3][3] = 1.0;
 
269
 
 
270
  matrix4_multiply(&tmp, matrix, matrix);
 
271
}
 
272
 
 
273
static void
 
274
matrix4_zrotate(RS_MATRIX4 *matrix, double rs, double rc)
 
275
{
 
276
  RS_MATRIX4 tmp;
 
277
 
 
278
  tmp.coeff[0][0] = rc;
 
279
  tmp.coeff[1][0] = rs;
 
280
  tmp.coeff[2][0] = 0.0;
 
281
  tmp.coeff[3][0] = 0.0;
 
282
 
 
283
  tmp.coeff[0][1] = -rs;
 
284
  tmp.coeff[1][1] = rc;
 
285
  tmp.coeff[2][1] = 0.0;
 
286
  tmp.coeff[3][1] = 0.0;
 
287
 
 
288
  tmp.coeff[0][2] = 0.0;
 
289
  tmp.coeff[1][2] = 0.0;
 
290
  tmp.coeff[2][2] = 1.0;
 
291
  tmp.coeff[3][2] = 0.0;
 
292
 
 
293
  tmp.coeff[0][3] = 0.0;
 
294
  tmp.coeff[1][3] = 0.0;
 
295
  tmp.coeff[2][3] = 0.0;
 
296
  tmp.coeff[3][3] = 1.0;
 
297
 
 
298
  matrix4_multiply(&tmp, matrix, matrix);
 
299
}
 
300
 
 
301
void
 
302
matrix4_color_normalize(RS_MATRIX4 *mat)
 
303
{
 
304
        int row,col;
 
305
        double lum;
 
306
 
 
307
        for(row=0;row<3;row++)
 
308
        {
 
309
                lum = 0.0;
 
310
                for(col=0;col<3;col++)
 
311
                        lum += mat->coeff[row][col];
 
312
                for(col=0;col<3;col++)
 
313
                        mat->coeff[row][col] /= lum;
 
314
        }
 
315
        return;
 
316
}
 
317
 
 
318
void
 
319
matrix4_color_saturate(RS_MATRIX4 *mat, double sat)
 
320
{
 
321
        RS_MATRIX4 tmp;
 
322
 
 
323
        if (sat == 1.0) return;
 
324
 
 
325
        tmp.coeff[0][0] = (1.0-sat)*RLUM + sat;
 
326
        tmp.coeff[1][0] = (1.0-sat)*RLUM;
 
327
        tmp.coeff[2][0] = (1.0-sat)*RLUM;
 
328
        tmp.coeff[3][0] = 0.0;
 
329
 
 
330
        tmp.coeff[0][1] = (1.0-sat)*GLUM;
 
331
        tmp.coeff[1][1] = (1.0-sat)*GLUM + sat;
 
332
        tmp.coeff[2][1] = (1.0-sat)*GLUM;
 
333
        tmp.coeff[3][1] = 0.0;
 
334
 
 
335
        tmp.coeff[0][2] = (1.0-sat)*BLUM;
 
336
        tmp.coeff[1][2] = (1.0-sat)*BLUM;
 
337
        tmp.coeff[2][2] = (1.0-sat)*BLUM + sat;
 
338
        tmp.coeff[3][2] = 0.0;
 
339
 
 
340
        tmp.coeff[0][3] = 0.0;
 
341
        tmp.coeff[1][3] = 0.0;
 
342
        tmp.coeff[2][3] = 0.0;
 
343
        tmp.coeff[3][3] = 1.0;
 
344
        matrix4_multiply(mat, &tmp, mat);
 
345
}
 
346
 
 
347
static void
 
348
matrix4_affine_transform_3dpoint(RS_MATRIX4 *matrix, double x, double y, double z, double *tx, double *ty, double *tz)
 
349
{
 
350
        *tx = x*matrix->coeff[0][0] + y*matrix->coeff[0][1]
 
351
                + z*matrix->coeff[0][2] + matrix->coeff[0][3];
 
352
        *ty = x*matrix->coeff[1][0] + y*matrix->coeff[1][1]
 
353
                + z*matrix->coeff[1][2] + matrix->coeff[1][3];
 
354
        *tz = x*matrix->coeff[2][0] + y*matrix->coeff[2][1]
 
355
                + z*matrix->coeff[2][2] + matrix->coeff[2][3];
 
356
}
 
357
 
 
358
void
 
359
matrix4_color_hue(RS_MATRIX4 *mat, double rot)
 
360
{
 
361
        RS_MATRIX4 tmp;
 
362
        double mag;
 
363
        double lx, ly, lz;
 
364
        double xrs, xrc;
 
365
        double yrs, yrc;
 
366
        double zrs, zrc;
 
367
        double zsx, zsy;
 
368
 
 
369
        if (rot==0.0) return;
 
370
 
 
371
        matrix4_identity(&tmp);
 
372
 
 
373
        /* rotate the grey vector into positive Z */
 
374
        mag = sqrt(2.0);
 
375
        xrs = 1.0/mag;
 
376
        xrc = 1.0/mag;
 
377
        matrix4_xrotate(&tmp, xrs, xrc);
 
378
 
 
379
        mag = sqrt(3.0);
 
380
        yrs = -1.0/mag;
 
381
        yrc = sqrt(2.0)/mag;
 
382
        matrix4_yrotate(&tmp, yrs ,yrc);
 
383
 
 
384
        /* shear the space to make the luminance plane horizontal */
 
385
        matrix4_affine_transform_3dpoint(&tmp,RLUM,GLUM,BLUM,&lx,&ly,&lz);
 
386
        zsx = lx/lz;
 
387
        zsy = ly/lz;
 
388
        matrix4_zshear(&tmp, zsx, zsy);
 
389
 
 
390
        /* rotate the hue */
 
391
        zrs = sin(rot*M_PI/180.0);
 
392
        zrc = cos(rot*M_PI/180.0);
 
393
        matrix4_zrotate(&tmp, zrs, zrc);
 
394
 
 
395
        /* unshear the space to put the luminance plane back */
 
396
        matrix4_zshear(&tmp, -zsx, -zsy);
 
397
 
 
398
        /* rotate the grey vector back into place */
 
399
        matrix4_yrotate(&tmp,-yrs,yrc);
 
400
        matrix4_xrotate(&tmp,-xrs,xrc);
 
401
        matrix4_multiply(mat,&tmp,mat);
 
402
}
 
403
 
 
404
void
 
405
matrix4_color_exposure(RS_MATRIX4 *mat, double exp)
 
406
{
 
407
        double expcom = pow(2.0, exp);
 
408
        mat->coeff[0][0] *= expcom;
 
409
        mat->coeff[0][1] *= expcom;
 
410
        mat->coeff[0][2] *= expcom;
 
411
        mat->coeff[1][0] *= expcom;
 
412
        mat->coeff[1][1] *= expcom;
 
413
        mat->coeff[1][2] *= expcom;
 
414
        mat->coeff[2][0] *= expcom;
 
415
        mat->coeff[2][1] *= expcom;
 
416
        mat->coeff[2][2] *= expcom;
 
417
        return;
 
418
}
 
419
 
 
420
void
 
421
matrix3_identity (RS_MATRIX3 *matrix)
 
422
{
 
423
  static const RS_MATRIX3 identity = { { { 1.0, 0.0, 0.0 },
 
424
                                          { 0.0, 1.0, 0.0 },
 
425
                                          { 0.0, 0.0, 1.0 } } };
 
426
 
 
427
  *matrix = identity;
 
428
}
 
429
 
 
430
RS_MATRIX3
 
431
matrix3_invert(const RS_MATRIX3 *matrix)
 
432
{
 
433
        int j,k;
 
434
 
 
435
        double a00 = matrix->coeff[0][0];
 
436
        double a01 = matrix->coeff[0][1];
 
437
        double a02 = matrix->coeff[0][2];
 
438
        double a10 = matrix->coeff[1][0];
 
439
        double a11 = matrix->coeff[1][1];
 
440
        double a12 = matrix->coeff[1][2];
 
441
        double a20 = matrix->coeff[2][0];
 
442
        double a21 = matrix->coeff[2][1];
 
443
        double a22 = matrix->coeff[2][2];
 
444
 
 
445
        double temp[3][3];
 
446
 
 
447
        temp[0][0] = a11 * a22 - a21 * a12;
 
448
        temp[0][1] = a21 * a02 - a01 * a22;
 
449
        temp[0][2] = a01 * a12 - a11 * a02;
 
450
        temp[1][0] = a20 * a12 - a10 * a22;
 
451
        temp[1][1] = a00 * a22 - a20 * a02;
 
452
        temp[1][2] = a10 * a02 - a00 * a12;
 
453
        temp[2][0] = a10 * a21 - a20 * a11;
 
454
        temp[2][1] = a20 * a01 - a00 * a21;
 
455
        temp[2][2] = a00 * a11 - a10 * a01;
 
456
 
 
457
        double det = (a00 * temp[0][0] + a01 * temp[1][0] + a02 * temp[2][0]);
 
458
 
 
459
        RS_MATRIX3 B;
 
460
 
 
461
        for (j = 0; j < 3; j++)
 
462
                for (k = 0; k < 3; k++)
 
463
                        B.coeff[j][k] = temp[j][k] / det;
 
464
 
 
465
        return B;
 
466
}
 
467
 
 
468
void
 
469
matrix3_multiply(const RS_MATRIX3 *left, const RS_MATRIX3 *right, RS_MATRIX3 *result)
 
470
{
 
471
  int i, j;
 
472
  RS_MATRIX3 tmp;
 
473
  double t1, t2, t3;
 
474
 
 
475
  for (i = 0; i < 3; i++)
 
476
    {
 
477
      t1 = left->coeff[i][0];
 
478
      t2 = left->coeff[i][1];
 
479
      t3 = left->coeff[i][2];
 
480
 
 
481
      for (j = 0; j < 3; j++)
 
482
        {
 
483
          tmp.coeff[i][j]  = t1 * right->coeff[0][j];
 
484
          tmp.coeff[i][j] += t2 * right->coeff[1][j];
 
485
          tmp.coeff[i][j] += t3 * right->coeff[2][j];
 
486
        }
 
487
    }
 
488
        *result = tmp;
 
489
        return;
 
490
}
 
491
 
 
492
RS_VECTOR3
 
493
vector3_multiply_matrix(const RS_VECTOR3 *vec, const RS_MATRIX3 *matrix)
 
494
{
 
495
        RS_VECTOR3 result;
 
496
 
 
497
        result.x = vec->x * matrix->coeff[0][0] + vec->y * matrix->coeff[0][1] + vec->z * matrix->coeff[0][2];
 
498
        result.y = vec->x * matrix->coeff[1][0] + vec->y * matrix->coeff[1][1] + vec->z * matrix->coeff[1][2];
 
499
        result.z = vec->x * matrix->coeff[2][0] + vec->y * matrix->coeff[2][1] + vec->z * matrix->coeff[2][2];
 
500
 
 
501
        return result;
 
502
}
 
503
 
 
504
float
 
505
matrix3_max(const RS_MATRIX3 *matrix)
 
506
{
 
507
        gfloat max = 0.0;
 
508
        int i, j;
 
509
        for (i = 0; i < 3; i++)
 
510
                for (j = 0; j < 3; j++)
 
511
                        max = MAX(max, matrix->coeff[i][j]);
 
512
 
 
513
        return max;
 
514
}
 
515
 
 
516
void
 
517
matrix3_scale(RS_MATRIX3 *matrix, const float scale, RS_MATRIX3 *result)
 
518
{
 
519
        int i, j;
 
520
        for (i = 0; i < 3; i++)
 
521
                for (j = 0; j < 3; j++)
 
522
                        result->coeff[i][j] = matrix->coeff[i][j] * scale;
 
523
}
 
524
 
 
525
void
 
526
matrix3_interpolate(const RS_MATRIX3 *a, const RS_MATRIX3 *b, const float alpha, RS_MATRIX3 *result)
 
527
{
 
528
#define LERP(a,b,g) (((b) - (a)) * (g) + (a))
 
529
        int i, j;
 
530
 
 
531
        for (i = 0; i < 3; i++)
 
532
                for (j = 0; j < 3; j++)
 
533
                        result->coeff[i][j] = LERP(a->coeff[i][j], b->coeff[i][j], alpha);
 
534
#undef LERP
 
535
}
 
536
 
 
537
float
 
538
matrix3_weight(const RS_MATRIX3 *mat)
 
539
{
 
540
        float weight = 
 
541
                mat->coeff[0][0]
 
542
                + mat->coeff[0][1]
 
543
                + mat->coeff[0][2]
 
544
                + mat->coeff[1][0]
 
545
                + mat->coeff[1][1]
 
546
                + mat->coeff[1][2]
 
547
                + mat->coeff[2][0]
 
548
                + mat->coeff[2][1]
 
549
                + mat->coeff[2][2];
 
550
        return(weight);
 
551
}
 
552
 
 
553
void
 
554
matrix3_to_matrix3int(RS_MATRIX3 *matrix, RS_MATRIX3Int *matrixi)
 
555
{
 
556
        int a,b;
 
557
        for(a=0;a<3;a++)
 
558
                for(b=0;b<3;b++)
 
559
                {
 
560
                        /* Check that a (unsigned 16 bit) x (matrix coeff shifted up RESOLUTION) fits within signed 32 bit value */
 
561
                        /* Adjust this if MATRIX_RESOLUTION is adjusted */
 
562
                        g_assert((matrix->coeff[a][b] < 16.0) && (matrix->coeff[a][b] > -16.0));
 
563
 
 
564
                        matrixi->coeff[a][b] = (int) (matrix->coeff[a][b] * (double) (1<<MATRIX_RESOLUTION));
 
565
                }
 
566
}
 
567
 
 
568
/*
 
569
 
 
570
Affine transformations
 
571
 
 
572
The following transformation matrix is used:
 
573
 
 
574
 [  a  b u ]
 
575
 [  c  d v ]
 
576
 [ tx ty w ]
 
577
 
 
578
a: x scale
 
579
b: y skew
 
580
c: x skew
 
581
d: y scale
 
582
tx: x translation (offset)
 
583
ty: y translation (offset)
 
584
 
 
585
(u, v & w unused)
 
586
 
 
587
x' = x*a + y*c + tx
 
588
y' = x*b + y*d + ty
 
589
 
 
590
*/
 
591
 
 
592
void
 
593
matrix3_affine_invert(RS_MATRIX3 *mat)
 
594
{
 
595
        RS_MATRIX3 tmp;
 
596
        const double reverse_det = 1.0/(mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
 
597
        matrix3_identity(&tmp);
 
598
        tmp.coeff[0][0] =  mat->coeff[1][1] * reverse_det;
 
599
        tmp.coeff[0][1] = -mat->coeff[0][1] * reverse_det;
 
600
        tmp.coeff[1][0] = -mat->coeff[1][0] * reverse_det;
 
601
        tmp.coeff[1][1] =  mat->coeff[0][0] * reverse_det;
 
602
        tmp.coeff[2][0] =  (mat->coeff[1][0]*mat->coeff[2][1] - mat->coeff[1][1]*mat->coeff[2][0])/
 
603
                (mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
 
604
        tmp.coeff[2][1] = -(mat->coeff[0][0]*mat->coeff[2][1] - mat->coeff[0][1]*mat->coeff[2][0])/
 
605
                (mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
 
606
        *mat = tmp;
 
607
}
 
608
 
 
609
void
 
610
matrix3_affine_scale(RS_MATRIX3 *matrix, double xscale, double yscale)
 
611
{
 
612
        RS_MATRIX3 tmp;
 
613
        matrix3_identity(&tmp);
 
614
        tmp.coeff[0][0] *= xscale;
 
615
        tmp.coeff[1][1] *= yscale;
 
616
        matrix3_multiply(matrix, &tmp, matrix);
 
617
        return;
 
618
}
 
619
 
 
620
void
 
621
matrix3_affine_translate(RS_MATRIX3 *matrix, double xtrans, double ytrans)
 
622
{
 
623
        matrix->coeff[2][0] += xtrans;
 
624
        matrix->coeff[2][1] += ytrans;
 
625
        return;
 
626
}
 
627
 
 
628
void
 
629
matrix3_affine_rotate(RS_MATRIX3 *matrix, double degrees)
 
630
{
 
631
        RS_MATRIX3 tmp;
 
632
        const double s = sin (degrees * M_PI / 180.0);
 
633
        const double c = cos (degrees * M_PI / 180.0);
 
634
 
 
635
        matrix3_identity(&tmp);
 
636
        tmp.coeff[0][0] = c;
 
637
        tmp.coeff[0][1] = s;
 
638
        tmp.coeff[1][0] = -s;
 
639
        tmp.coeff[1][1] = c;
 
640
        matrix3_multiply(matrix, &tmp, matrix);
 
641
        return;
 
642
}
 
643
 
 
644
inline void
 
645
matrix3_affine_transform_point(RS_MATRIX3 *matrix, double x, double y, double *x2, double *y2)
 
646
{
 
647
        const double x_tmp = x*matrix->coeff[0][0] + y*matrix->coeff[1][0] + matrix->coeff[2][0];
 
648
        const double y_tmp = x*matrix->coeff[0][1] + y*matrix->coeff[1][1] + matrix->coeff[2][1];
 
649
        *x2 = x_tmp;
 
650
        *y2 = y_tmp;
 
651
        return;
 
652
}
 
653
 
 
654
inline void
 
655
matrix3_affine_transform_point_int(RS_MATRIX3 *matrix, int x, int y, int *x2, int *y2)
 
656
{
 
657
        const int x_tmp = (int) (x*matrix->coeff[0][0] + y*matrix->coeff[1][0] + matrix->coeff[2][0]);
 
658
        const int y_tmp = (int) (x*matrix->coeff[0][1] + y*matrix->coeff[1][1] + matrix->coeff[2][1]);
 
659
        *x2 = x_tmp;
 
660
        *y2 = y_tmp;
 
661
        return;
 
662
}
 
663
 
 
664
void
 
665
matrix3_affine_get_minmax(RS_MATRIX3 *matrix,
 
666
        double *minx, double *miny,
 
667
        double *maxx, double *maxy,
 
668
        double x1, double y1,
 
669
        double x2, double y2)
 
670
{
 
671
        double x,y;
 
672
        *minx = *miny = 100000.0;
 
673
        *maxx = *maxy = 0.0;
 
674
 
 
675
        matrix3_affine_transform_point(matrix, x1, y1, &x, &y);
 
676
        if (x<*minx) *minx=x;
 
677
        if (x>*maxx) *maxx=x;
 
678
        if (y<*miny) *miny=y;
 
679
        if (y>*maxy) *maxy=y;
 
680
        matrix3_affine_transform_point(matrix, x2, y1, &x, &y);
 
681
        if (x<*minx) *minx=x;
 
682
        if (x>*maxx) *maxx=x;
 
683
        if (y<*miny) *miny=y;
 
684
        if (y>*maxy) *maxy=y;
 
685
        matrix3_affine_transform_point(matrix, x1, y2, &x, &y);
 
686
        if (x<*minx) *minx=x;
 
687
        if (x>*maxx) *maxx=x;
 
688
        if (y<*miny) *miny=y;
 
689
        if (y>*maxy) *maxy=y;
 
690
        matrix3_affine_transform_point(matrix, x2, y2, &x, &y);
 
691
        if (x<*minx) *minx=x;
 
692
        if (x>*maxx) *maxx=x;
 
693
        if (y<*miny) *miny=y;
 
694
        if (y>*maxy) *maxy=y;
 
695
 
 
696
        return;
 
697
}
 
698
 
 
699
/**
 
700
 * Interpolate an array of unsigned integers
 
701
 * @param input_dataset An array of unsigned integers to be inperpolated
 
702
 * @param input_length The length of the input array
 
703
 * @param output_dataset An array of unsigned integers for output or NULL
 
704
 * @param output_length The length of the output array
 
705
 * @param max A pointer to an unsigned int or NULL
 
706
 * @return the interpolated dataset
 
707
 */
 
708
unsigned int *
 
709
interpolate_dataset_int(unsigned int *input_dataset, unsigned int input_length, unsigned int *output_dataset, unsigned int output_length, unsigned int *max)
 
710
{
 
711
        const double scale = ((double)input_length) / ((double)output_length);
 
712
        int i, source1, source2;
 
713
        float source;
 
714
        float weight1, weight2;
 
715
 
 
716
        if (output_dataset == NULL)
 
717
                output_dataset = malloc(sizeof(unsigned int)*output_length);
 
718
 
 
719
        for(i=0;i<output_length;i++)
 
720
        {
 
721
                source = ((float)i) * scale;
 
722
                source1 = (int) floor(source);
 
723
                source2 = source1+1;
 
724
                weight1 = 1.0f - (source - floor(source));
 
725
                weight2 = 1.0f - weight1;
 
726
                output_dataset[i] = (unsigned int) (((float)input_dataset[source1]) * weight1
 
727
                                        + ((float)input_dataset[source2]) * weight2);
 
728
                if (max)
 
729
                        if (output_dataset[i] > (*max))
 
730
                                *max = output_dataset[i];
 
731
        }
 
732
 
 
733
        return output_dataset;
 
734
}