← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/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:
src/rs-math.c
1
 
/*
2
 
 * Copyright (C) 2006-2009 Anders Brander <anders@brander.dk> and 
3
 
 * Anders Kvist <akv@lnxbx.dk>
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 "color.h"
23
 
#include "rs-math.h"
24
 
 
25
 
static void matrix4_zshear (RS_MATRIX4 *matrix, double dx, double dy);
26
 
static void matrix4_xrotate(RS_MATRIX4 *matrix, double rs, double rc);
27
 
static void matrix4_yrotate(RS_MATRIX4 *matrix, double rs, double rc);
28
 
static void matrix4_zrotate(RS_MATRIX4 *matrix, double rs, double rc);
29
 
static void matrix4_affine_transform_3dpoint(RS_MATRIX4 *matrix, double x, double y, double z, double *tx, double *ty, double *tz);
30
 
 
31
 
void
32
 
printmat(RS_MATRIX4 *mat)
33
 
{
34
 
        int x, y;
35
 
 
36
 
        for(y=0; y<4; y++)
37
 
        {
38
 
                for(x=0; x<4; x++)
39
 
                        printf("%f ",mat->coeff[x][y]);
40
 
                printf("\n");
41
 
        }
42
 
        printf("\n");
43
 
}
44
 
 
45
 
void
46
 
matrix4_identity (RS_MATRIX4 *matrix)
47
 
{
48
 
  static const RS_MATRIX4 identity = { { { 1.0, 0.0, 0.0, 0.0 },
49
 
                                          { 0.0, 1.0, 0.0, 0.0 },
50
 
                                          { 0.0, 0.0, 1.0, 0.0 },
51
 
                                          { 0.0, 0.0, 0.0, 1.0 } } };
52
 
 
53
 
  *matrix = identity;
54
 
}
55
 
 
56
 
void
57
 
matrix4_multiply(const RS_MATRIX4 *left, RS_MATRIX4 *right, RS_MATRIX4 *result)
58
 
{
59
 
  int i, j;
60
 
  RS_MATRIX4 tmp;
61
 
  double t1, t2, t3, t4;
62
 
 
63
 
  for (i = 0; i < 4; i++)
64
 
    {
65
 
      t1 = left->coeff[i][0];
66
 
      t2 = left->coeff[i][1];
67
 
      t3 = left->coeff[i][2];
68
 
      t4 = left->coeff[i][3];
69
 
 
70
 
      for (j = 0; j < 4; j++)
71
 
        {
72
 
          tmp.coeff[i][j]  = t1 * right->coeff[0][j];
73
 
          tmp.coeff[i][j] += t2 * right->coeff[1][j];
74
 
          tmp.coeff[i][j] += t3 * right->coeff[2][j];
75
 
          tmp.coeff[i][j] += t4 * right->coeff[3][j];
76
 
        }
77
 
    }
78
 
  *result = tmp;
79
 
}
80
 
 
81
 
/* copied almost verbatim from  dcraw.c:pseudoinverse()
82
 
 - but this one doesn't transpose */
83
 
void
84
 
matrix4_color_invert(const RS_MATRIX4 *in, RS_MATRIX4 *out)
85
 
{
86
 
        RS_MATRIX4 tmp;
87
 
        double work[3][6], num;
88
 
        int i, j, k;
89
 
 
90
 
        matrix4_identity(&tmp);
91
 
 
92
 
        for (i=0; i < 3; i++)
93
 
        {
94
 
                for (j=0; j < 6; j++)
95
 
                        work[i][j] = j == i+3;
96
 
                for (j=0; j < 3; j++)
97
 
                        for (k=0; k < 3; k++)
98
 
                                work[i][j] += in->coeff[k][i] * in->coeff[k][j];
99
 
        }
100
 
        for (i=0; i < 3; i++)
101
 
        {
102
 
                num = work[i][i];
103
 
                for (j=0; j < 6; j++)
104
 
                        work[i][j] /= num;
105
 
                for (k=0; k < 3; k++)
106
 
                {
107
 
                        if (k==i)
108
 
                                continue;
109
 
                        num = work[k][i];
110
 
                        for (j=0; j < 6; j++)
111
 
                                work[k][j] -= work[i][j] * num;
112
 
                }
113
 
        }
114
 
        for (i=0; i < 3; i++)
115
 
                for (j=0; j < 3; j++)
116
 
                        for (tmp.coeff[i][j]=k=0; k < 3; k++)
117
 
                                tmp.coeff[i][j] += work[j][k+3] * in->coeff[i][k];
118
 
        for (i=0; i < 4; i++)
119
 
                for (j=0; j < 4; j++)
120
 
                        out->coeff[i][j] = tmp.coeff[j][i];
121
 
}
122
 
 
123
 
static void
124
 
matrix4_zshear (RS_MATRIX4 *matrix, double dx, double dy)
125
 
{
126
 
  RS_MATRIX4 zshear;
127
 
 
128
 
  zshear.coeff[0][0] = 1.0;
129
 
  zshear.coeff[1][0] = 0.0;
130
 
  zshear.coeff[2][0] = dx;
131
 
  zshear.coeff[3][0] = 0.0;
132
 
 
133
 
  zshear.coeff[0][1] = 0.0;
134
 
  zshear.coeff[1][1] = 1.0;
135
 
  zshear.coeff[2][1] = dy;
136
 
  zshear.coeff[3][1] = 0.0;
137
 
 
138
 
  zshear.coeff[0][2] = 0.0;
139
 
  zshear.coeff[1][2] = 0.0;
140
 
  zshear.coeff[2][2] = 1.0;
141
 
  zshear.coeff[3][2] = 0.0;
142
 
 
143
 
  zshear.coeff[0][3] = 0.0;
144
 
  zshear.coeff[1][3] = 0.0;
145
 
  zshear.coeff[2][3] = 0.0;
146
 
  zshear.coeff[3][3] = 1.0;
147
 
 
148
 
  matrix4_multiply(&zshear, matrix, matrix);
149
 
}
150
 
 
151
 
void
152
 
matrix4_to_matrix4int(RS_MATRIX4 *matrix, RS_MATRIX4Int *matrixi)
153
 
{
154
 
  int a,b;
155
 
  for(a=0;a<4;a++)
156
 
    for(b=0;b<4;b++)
157
 
      matrixi->coeff[a][b] = (int) (matrix->coeff[a][b] * (double) (1<<MATRIX_RESOLUTION));
158
 
  return;
159
 
}
160
 
 
161
 
static void
162
 
matrix4_xrotate(RS_MATRIX4 *matrix, double rs, double rc)
163
 
{
164
 
  RS_MATRIX4 tmp;
165
 
 
166
 
  tmp.coeff[0][0] = 1.0;
167
 
  tmp.coeff[1][0] = 0.0;
168
 
  tmp.coeff[2][0] = 0.0;
169
 
  tmp.coeff[3][0] = 0.0;
170
 
 
171
 
  tmp.coeff[0][1] = 0.0;
172
 
  tmp.coeff[1][1] = rc;
173
 
  tmp.coeff[2][1] = rs;
174
 
  tmp.coeff[3][1] = 0.0;
175
 
 
176
 
  tmp.coeff[0][2] = 0.0;
177
 
  tmp.coeff[1][2] = -rs;
178
 
  tmp.coeff[2][2] = rc;
179
 
  tmp.coeff[3][2] = 0.0;
180
 
 
181
 
  tmp.coeff[0][3] = 0.0;
182
 
  tmp.coeff[1][3] = 0.0;
183
 
  tmp.coeff[2][3] = 0.0;
184
 
  tmp.coeff[3][3] = 1.0;
185
 
 
186
 
  matrix4_multiply(&tmp, matrix, matrix);
187
 
}
188
 
 
189
 
static void
190
 
matrix4_yrotate(RS_MATRIX4 *matrix, double rs, double rc)
191
 
{
192
 
  RS_MATRIX4 tmp;
193
 
 
194
 
  tmp.coeff[0][0] = rc;
195
 
  tmp.coeff[1][0] = 0.0;
196
 
  tmp.coeff[2][0] = -rs;
197
 
  tmp.coeff[3][0] = 0.0;
198
 
 
199
 
  tmp.coeff[0][1] = 0.0;
200
 
  tmp.coeff[1][1] = 1.0;
201
 
  tmp.coeff[2][1] = 0.0;
202
 
  tmp.coeff[3][1] = 0.0;
203
 
 
204
 
  tmp.coeff[0][2] = rs;
205
 
  tmp.coeff[1][2] = 0.0;
206
 
  tmp.coeff[2][2] = rc;
207
 
  tmp.coeff[3][2] = 0.0;
208
 
 
209
 
  tmp.coeff[0][3] = 0.0;
210
 
  tmp.coeff[1][3] = 0.0;
211
 
  tmp.coeff[2][3] = 0.0;
212
 
  tmp.coeff[3][3] = 1.0;
213
 
 
214
 
  matrix4_multiply(&tmp, matrix, matrix);
215
 
}
216
 
 
217
 
static void
218
 
matrix4_zrotate(RS_MATRIX4 *matrix, double rs, double rc)
219
 
{
220
 
  RS_MATRIX4 tmp;
221
 
 
222
 
  tmp.coeff[0][0] = rc;
223
 
  tmp.coeff[1][0] = rs;
224
 
  tmp.coeff[2][0] = 0.0;
225
 
  tmp.coeff[3][0] = 0.0;
226
 
 
227
 
  tmp.coeff[0][1] = -rs;
228
 
  tmp.coeff[1][1] = rc;
229
 
  tmp.coeff[2][1] = 0.0;
230
 
  tmp.coeff[3][1] = 0.0;
231
 
 
232
 
  tmp.coeff[0][2] = 0.0;
233
 
  tmp.coeff[1][2] = 0.0;
234
 
  tmp.coeff[2][2] = 1.0;
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
 
void
246
 
matrix4_color_normalize(RS_MATRIX4 *mat)
247
 
{
248
 
        int row,col;
249
 
        double lum;
250
 
 
251
 
        for(row=0;row<3;row++)
252
 
        {
253
 
                lum = 0.0;
254
 
                for(col=0;col<3;col++)
255
 
                        lum += mat->coeff[row][col];
256
 
                for(col=0;col<3;col++)
257
 
                        mat->coeff[row][col] /= lum;
258
 
        }
259
 
        return;
260
 
}
261
 
 
262
 
void
263
 
matrix4_color_saturate(RS_MATRIX4 *mat, double sat)
264
 
{
265
 
        RS_MATRIX4 tmp;
266
 
 
267
 
        if (sat == 1.0) return;
268
 
 
269
 
        tmp.coeff[0][0] = (1.0-sat)*RLUM + sat;
270
 
        tmp.coeff[1][0] = (1.0-sat)*RLUM;
271
 
        tmp.coeff[2][0] = (1.0-sat)*RLUM;
272
 
        tmp.coeff[3][0] = 0.0;
273
 
 
274
 
        tmp.coeff[0][1] = (1.0-sat)*GLUM;
275
 
        tmp.coeff[1][1] = (1.0-sat)*GLUM + sat;
276
 
        tmp.coeff[2][1] = (1.0-sat)*GLUM;
277
 
        tmp.coeff[3][1] = 0.0;
278
 
 
279
 
        tmp.coeff[0][2] = (1.0-sat)*BLUM;
280
 
        tmp.coeff[1][2] = (1.0-sat)*BLUM;
281
 
        tmp.coeff[2][2] = (1.0-sat)*BLUM + sat;
282
 
        tmp.coeff[3][2] = 0.0;
283
 
 
284
 
        tmp.coeff[0][3] = 0.0;
285
 
        tmp.coeff[1][3] = 0.0;
286
 
        tmp.coeff[2][3] = 0.0;
287
 
        tmp.coeff[3][3] = 1.0;
288
 
        matrix4_multiply(mat, &tmp, mat);
289
 
}
290
 
 
291
 
static void
292
 
matrix4_affine_transform_3dpoint(RS_MATRIX4 *matrix, double x, double y, double z, double *tx, double *ty, double *tz)
293
 
{
294
 
        *tx = x*matrix->coeff[0][0] + y*matrix->coeff[0][1]
295
 
                + z*matrix->coeff[0][2] + matrix->coeff[0][3];
296
 
        *ty = x*matrix->coeff[1][0] + y*matrix->coeff[1][1]
297
 
                + z*matrix->coeff[1][2] + matrix->coeff[1][3];
298
 
        *tz = x*matrix->coeff[2][0] + y*matrix->coeff[2][1]
299
 
                + z*matrix->coeff[2][2] + matrix->coeff[2][3];
300
 
}
301
 
 
302
 
void
303
 
matrix4_color_hue(RS_MATRIX4 *mat, double rot)
304
 
{
305
 
        RS_MATRIX4 tmp;
306
 
        double mag;
307
 
        double lx, ly, lz;
308
 
        double xrs, xrc;
309
 
        double yrs, yrc;
310
 
        double zrs, zrc;
311
 
        double zsx, zsy;
312
 
 
313
 
        if (rot==0.0) return;
314
 
 
315
 
        matrix4_identity(&tmp);
316
 
 
317
 
        /* rotate the grey vector into positive Z */
318
 
        mag = sqrt(2.0);
319
 
        xrs = 1.0/mag;
320
 
        xrc = 1.0/mag;
321
 
        matrix4_xrotate(&tmp, xrs, xrc);
322
 
 
323
 
        mag = sqrt(3.0);
324
 
        yrs = -1.0/mag;
325
 
        yrc = sqrt(2.0)/mag;
326
 
        matrix4_yrotate(&tmp, yrs ,yrc);
327
 
 
328
 
        /* shear the space to make the luminance plane horizontal */
329
 
        matrix4_affine_transform_3dpoint(&tmp,RLUM,GLUM,BLUM,&lx,&ly,&lz);
330
 
        zsx = lx/lz;
331
 
        zsy = ly/lz;
332
 
        matrix4_zshear(&tmp, zsx, zsy);
333
 
 
334
 
        /* rotate the hue */
335
 
        zrs = sin(rot*M_PI/180.0);
336
 
        zrc = cos(rot*M_PI/180.0);
337
 
        matrix4_zrotate(&tmp, zrs, zrc);
338
 
 
339
 
        /* unshear the space to put the luminance plane back */
340
 
        matrix4_zshear(&tmp, -zsx, -zsy);
341
 
 
342
 
        /* rotate the grey vector back into place */
343
 
        matrix4_yrotate(&tmp,-yrs,yrc);
344
 
        matrix4_xrotate(&tmp,-xrs,xrc);
345
 
        matrix4_multiply(mat,&tmp,mat);
346
 
}
347
 
 
348
 
void
349
 
matrix4_color_exposure(RS_MATRIX4 *mat, double exp)
350
 
{
351
 
        double expcom = pow(2.0, exp);
352
 
        mat->coeff[0][0] *= expcom;
353
 
        mat->coeff[0][1] *= expcom;
354
 
        mat->coeff[0][2] *= expcom;
355
 
        mat->coeff[1][0] *= expcom;
356
 
        mat->coeff[1][1] *= expcom;
357
 
        mat->coeff[1][2] *= expcom;
358
 
        mat->coeff[2][0] *= expcom;
359
 
        mat->coeff[2][1] *= expcom;
360
 
        mat->coeff[2][2] *= expcom;
361
 
        return;
362
 
}
363
 
 
364
 
void
365
 
matrix3_identity (RS_MATRIX3 *matrix)
366
 
{
367
 
  static const RS_MATRIX3 identity = { { { 1.0, 0.0, 0.0 },
368
 
                                          { 0.0, 1.0, 0.0 },
369
 
                                          { 0.0, 0.0, 1.0 } } };
370
 
 
371
 
  *matrix = identity;
372
 
}
373
 
 
374
 
static void
375
 
matrix3_multiply(const RS_MATRIX3 *left, RS_MATRIX3 *right, RS_MATRIX3 *result)
376
 
{
377
 
  int i, j;
378
 
  RS_MATRIX3 tmp;
379
 
  double t1, t2, t3;
380
 
 
381
 
  for (i = 0; i < 3; i++)
382
 
    {
383
 
      t1 = left->coeff[i][0];
384
 
      t2 = left->coeff[i][1];
385
 
      t3 = left->coeff[i][2];
386
 
 
387
 
      for (j = 0; j < 3; j++)
388
 
        {
389
 
          tmp.coeff[i][j]  = t1 * right->coeff[0][j];
390
 
          tmp.coeff[i][j] += t2 * right->coeff[1][j];
391
 
          tmp.coeff[i][j] += t3 * right->coeff[2][j];
392
 
        }
393
 
    }
394
 
        *result = tmp;
395
 
        return;
396
 
}
397
 
 
398
 
float
399
 
matrix3_weight(const RS_MATRIX3 *mat)
400
 
{
401
 
        float weight = 
402
 
                mat->coeff[0][0]
403
 
                + mat->coeff[0][1]
404
 
                + mat->coeff[0][2]
405
 
                + mat->coeff[1][0]
406
 
                + mat->coeff[1][1]
407
 
                + mat->coeff[1][2]
408
 
                + mat->coeff[2][0]
409
 
                + mat->coeff[2][1]
410
 
                + mat->coeff[2][2];
411
 
        return(weight);
412
 
}
413
 
 
414
 
void
415
 
matrix3_to_matrix3int(RS_MATRIX3 *matrix, RS_MATRIX3Int *matrixi)
416
 
{
417
 
  int a,b;
418
 
  for(a=0;a<3;a++)
419
 
    for(b=0;b<3;b++)
420
 
      matrixi->coeff[a][b] = (int) (matrix->coeff[a][b] * (double) (1<<MATRIX_RESOLUTION));
421
 
  return;
422
 
}
423
 
 
424
 
/*
425
 
 
426
 
Affine transformations
427
 
 
428
 
The following transformation matrix is used:
429
 
 
430
 
 [  a  b u ]
431
 
 [  c  d v ]
432
 
 [ tx ty w ]
433
 
 
434
 
a: x scale
435
 
b: y skew
436
 
c: x skew
437
 
d: y scale
438
 
tx: x translation (offset)
439
 
ty: y translation (offset)
440
 
 
441
 
(u, v & w unused)
442
 
 
443
 
x' = x*a + y*c + tx
444
 
y' = x*b + y*d + ty
445
 
 
446
 
*/
447
 
 
448
 
void
449
 
matrix3_affine_invert(RS_MATRIX3 *mat)
450
 
{
451
 
        RS_MATRIX3 tmp;
452
 
        const double reverse_det = 1.0/(mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
453
 
        matrix3_identity(&tmp);
454
 
        tmp.coeff[0][0] =  mat->coeff[1][1] * reverse_det;
455
 
        tmp.coeff[0][1] = -mat->coeff[0][1] * reverse_det;
456
 
        tmp.coeff[1][0] = -mat->coeff[1][0] * reverse_det;
457
 
        tmp.coeff[1][1] =  mat->coeff[0][0] * reverse_det;
458
 
        tmp.coeff[2][0] =  (mat->coeff[1][0]*mat->coeff[2][1] - mat->coeff[1][1]*mat->coeff[2][0])/
459
 
                (mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
460
 
        tmp.coeff[2][1] = -(mat->coeff[0][0]*mat->coeff[2][1] - mat->coeff[0][1]*mat->coeff[2][0])/
461
 
                (mat->coeff[0][0]*mat->coeff[1][1] - mat->coeff[0][1]*mat->coeff[1][0]);
462
 
        *mat = tmp;
463
 
}
464
 
 
465
 
void
466
 
matrix3_affine_scale(RS_MATRIX3 *matrix, double xscale, double yscale)
467
 
{
468
 
        RS_MATRIX3 tmp;
469
 
        matrix3_identity(&tmp);
470
 
        tmp.coeff[0][0] *= xscale;
471
 
        tmp.coeff[1][1] *= yscale;
472
 
        matrix3_multiply(matrix, &tmp, matrix);
473
 
        return;
474
 
}
475
 
 
476
 
void
477
 
matrix3_affine_translate(RS_MATRIX3 *matrix, double xtrans, double ytrans)
478
 
{
479
 
        matrix->coeff[2][0] += xtrans;
480
 
        matrix->coeff[2][1] += ytrans;
481
 
        return;
482
 
}
483
 
 
484
 
void
485
 
matrix3_affine_rotate(RS_MATRIX3 *matrix, double degrees)
486
 
{
487
 
        RS_MATRIX3 tmp;
488
 
        const double s = sin (degrees * M_PI / 180.0);
489
 
        const double c = cos (degrees * M_PI / 180.0);
490
 
 
491
 
        matrix3_identity(&tmp);
492
 
        tmp.coeff[0][0] = c;
493
 
        tmp.coeff[0][1] = s;
494
 
        tmp.coeff[1][0] = -s;
495
 
        tmp.coeff[1][1] = c;
496
 
        matrix3_multiply(matrix, &tmp, matrix);
497
 
        return;
498
 
}
499
 
 
500
 
inline void
501
 
matrix3_affine_transform_point(RS_MATRIX3 *matrix, double x, double y, double *x2, double *y2)
502
 
{
503
 
        const double x_tmp = x*matrix->coeff[0][0] + y*matrix->coeff[1][0] + matrix->coeff[2][0];
504
 
        const double y_tmp = x*matrix->coeff[0][1] + y*matrix->coeff[1][1] + matrix->coeff[2][1];
505
 
        *x2 = x_tmp;
506
 
        *y2 = y_tmp;
507
 
        return;
508
 
}
509
 
 
510
 
inline void
511
 
matrix3_affine_transform_point_int(RS_MATRIX3 *matrix, int x, int y, int *x2, int *y2)
512
 
{
513
 
        const int x_tmp = (int) (x*matrix->coeff[0][0] + y*matrix->coeff[1][0] + matrix->coeff[2][0]);
514
 
        const int y_tmp = (int) (x*matrix->coeff[0][1] + y*matrix->coeff[1][1] + matrix->coeff[2][1]);
515
 
        *x2 = x_tmp;
516
 
        *y2 = y_tmp;
517
 
        return;
518
 
}
519
 
 
520
 
void
521
 
matrix3_affine_get_minmax(RS_MATRIX3 *matrix,
522
 
        double *minx, double *miny,
523
 
        double *maxx, double *maxy,
524
 
        double x1, double y1,
525
 
        double x2, double y2)
526
 
{
527
 
        double x,y;
528
 
        *minx = *miny = 100000.0;
529
 
        *maxx = *maxy = 0.0;
530
 
 
531
 
        matrix3_affine_transform_point(matrix, x1, y1, &x, &y);
532
 
        if (x<*minx) *minx=x;
533
 
        if (x>*maxx) *maxx=x;
534
 
        if (y<*miny) *miny=y;
535
 
        if (y>*maxy) *maxy=y;
536
 
        matrix3_affine_transform_point(matrix, x2, y1, &x, &y);
537
 
        if (x<*minx) *minx=x;
538
 
        if (x>*maxx) *maxx=x;
539
 
        if (y<*miny) *miny=y;
540
 
        if (y>*maxy) *maxy=y;
541
 
        matrix3_affine_transform_point(matrix, x1, y2, &x, &y);
542
 
        if (x<*minx) *minx=x;
543
 
        if (x>*maxx) *maxx=x;
544
 
        if (y<*miny) *miny=y;
545
 
        if (y>*maxy) *maxy=y;
546
 
        matrix3_affine_transform_point(matrix, x2, y2, &x, &y);
547
 
        if (x<*minx) *minx=x;
548
 
        if (x>*maxx) *maxx=x;
549
 
        if (y<*miny) *miny=y;
550
 
        if (y>*maxy) *maxy=y;
551
 
 
552
 
        return;
553
 
}
554
 
 
555
 
/**
556
 
 * Interpolate an array of unsigned integers
557
 
 * @param input_dataset An array of unsigned integers to be inperpolated
558
 
 * @param input_length The length of the input array
559
 
 * @param output_dataset An array of unsigned integers for output or NULL
560
 
 * @param output_length The length of the output array
561
 
 * @param max A pointer to an unsigned int or NULL
562
 
 * @return the interpolated dataset
563
 
 */
564
 
unsigned int *
565
 
interpolate_dataset_int(unsigned int *input_dataset, unsigned int input_length, unsigned int *output_dataset, unsigned int output_length, unsigned int *max)
566
 
{
567
 
        const double scale = ((double)input_length) / ((double)output_length);
568
 
        int i, source1, source2;
569
 
        float source;
570
 
        float weight1, weight2;
571
 
 
572
 
        if (output_dataset == NULL)
573
 
                output_dataset = malloc(sizeof(unsigned int)*output_length);
574
 
 
575
 
        for(i=0;i<output_length;i++)
576
 
        {
577
 
                source = ((float)i) * scale;
578
 
                source1 = (int) floor(source);
579
 
                source2 = source1+1;
580
 
                weight1 = 1.0f - (source - floor(source));
581
 
                weight2 = 1.0f - weight1;
582
 
                output_dataset[i] = (unsigned int) (((float)input_dataset[source1]) * weight1
583
 
                                        + ((float)input_dataset[source2]) * weight2);
584
 
                if (max)
585
 
                        if (output_dataset[i] > (*max))
586
 
                                *max = output_dataset[i];
587
 
        }
588
 
 
589
 
        return output_dataset;
590
 
}