← Back to branch summary

~ubuntu-branches/ubuntu/saucy/digikam/saucy

« back to all changes in this revision

Viewing changes to libs/dimg/filters/sharp/matrix.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christian Mangold
  • Date: 2010-04-09 21:30:01 UTC
  • mfrom: (1.2.28 upstream)
  • Revision ID: james.westby@ubuntu.com-20100409213001-4bfyibrd359rn7o3
Tags: 2:1.2.0-0ubuntu1
https://launchpad.net/bugs/560576
* New upstream release (LP: #560576)
* Remove all patches, fixed upstream
  - Remove quilt build-depend

Show diffs side-by-side

added added

removed removed

Lines of Context:
libs/dimg/filters/sharp/matrix.cpp
 
1
/* ============================================================
 
2
 *
 
3
 * This file is a part of digiKam project
 
4
 * http://www.digikam.org
 
5
 *
 
6
 * Date        : 2005-04-29
 
7
 * Description : refocus deconvolution matrix implementation.
 
8
 *
 
9
 * Copyright (C) 2005-2010 by Gilles Caulier <caulier dot gilles at gmail dot com>
 
10
 *
 
11
 * Original implementation from Refocus Gimp plug-in
 
12
 * Copyright (C) 1999-2003 Ernst Lippe
 
13
 *
 
14
 * This program is free software; you can redistribute it
 
15
 * and/or modify it under the terms of the GNU General
 
16
 * Public License as published by the Free Software Foundation;
 
17
 * either version 2, or (at your option)
 
18
 * any later version.
 
19
 *
 
20
 * This program is distributed in the hope that it will be useful,
 
21
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
22
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
23
 * GNU General Public License for more details.
 
24
 *
 
25
 * ============================================================ */
 
26
 
 
27
// Uncomment this line to debug matrix computation.
 
28
//#define RF_DEBUG 1
 
29
 
 
30
// Square
 
31
#define SQR(x) ((x) * (x))
 
32
 
 
33
#include "matrix.h"
 
34
 
 
35
// C ANSI includes
 
36
 
 
37
extern "C"
 
38
{
 
39
#include "f2c.h"
 
40
#include "clapack.h"
 
41
}
 
42
 
 
43
// C++ includes
 
44
 
 
45
#include <cmath>
 
46
 
 
47
// Qt includes
 
48
 
 
49
#include <QString>
 
50
 
 
51
// KDE includes
 
52
 
 
53
#include <kdebug.h>
 
54
 
 
55
namespace Digikam
 
56
{
 
57
 
 
58
Mat *RefocusMatrix::allocate_matrix (int nrows, int ncols)
 
59
{
 
60
    Mat *result = new Mat;
 
61
    memset (result, 0, sizeof(result));
 
62
 
 
63
    result->cols = ncols;
 
64
    result->rows = nrows;
 
65
    result->data = new double [nrows * ncols];
 
66
    memset (result->data, 0, nrows * ncols * sizeof(double));
 
67
 
 
68
    return (result);
 
69
}
 
70
 
 
71
void RefocusMatrix::finish_matrix (Mat * mat)
 
72
{
 
73
    delete [] mat->data;
 
74
}
 
75
 
 
76
void RefocusMatrix::finish_and_free_matrix (Mat * mat)
 
77
{
 
78
    delete [] mat->data;
 
79
    delete mat;
 
80
}
 
81
 
 
82
double *RefocusMatrix::mat_eltptr (Mat * mat, const int r, const int c)
 
83
{
 
84
    Q_ASSERT ((r >= 0) && (r < mat->rows));
 
85
    Q_ASSERT ((c >= 0) && (c < mat->rows));
 
86
    return (&(mat->data[mat->rows * c + r]));
 
87
}
 
88
 
 
89
double RefocusMatrix::mat_elt (const Mat * mat, const int r, const int c)
 
90
{
 
91
    Q_ASSERT ((r >= 0) && (r < mat->rows));
 
92
    Q_ASSERT ((c >= 0) && (c < mat->rows));
 
93
    return (mat->data[mat->rows * c + r]);
 
94
}
 
95
 
 
96
void RefocusMatrix::init_c_mat (CMat * mat, const int radius)
 
97
{
 
98
    mat->radius = radius;
 
99
    mat->row_stride = 2 * radius + 1;
 
100
    mat->data = new double [SQR (mat->row_stride)];
 
101
    memset (mat->data, 0, SQR (mat->row_stride) * sizeof(double));
 
102
    mat->center = mat->data + mat->row_stride * mat->radius + mat->radius;
 
103
}
 
104
 
 
105
CMat *RefocusMatrix::allocate_c_mat (const int radius)
 
106
{
 
107
    CMat *result = new CMat;
 
108
    memset(result, 0, sizeof(result));
 
109
    init_c_mat (result, radius);
 
110
    return (result);
 
111
}
 
112
 
 
113
void RefocusMatrix::finish_c_mat (CMat * mat)
 
114
{
 
115
    delete [] mat->data;
 
116
    mat->data = NULL;
 
117
}
 
118
 
 
119
inline double *RefocusMatrix::c_mat_eltptr (CMat * mat, const int col, const int row)
 
120
{
 
121
    Q_ASSERT ((qAbs (row) <= mat->radius) && (qAbs (col) <= mat->radius));
 
122
    return (mat->center + mat->row_stride * row + col);
 
123
}
 
124
 
 
125
inline double RefocusMatrix::c_mat_elt (const CMat * const mat, const int col, const int row)
 
126
{
 
127
    Q_ASSERT ((qAbs (row) <= mat->radius) && (qAbs (col) <= mat->radius));
 
128
    return (mat->center[mat->row_stride * row + col]);
 
129
}
 
130
 
 
131
void RefocusMatrix::convolve_mat (CMat * result, const CMat * const mata, const CMat * const matb)
 
132
{
 
133
    register int xr, yr, xa, ya;
 
134
 
 
135
    for (yr = -result->radius; yr <= result->radius; ++yr)
 
136
    {
 
137
        for (xr = -result->radius; xr <= result->radius; ++xr)
 
138
        {
 
139
            const int ya_low  = qMax (-mata->radius, yr - matb->radius);
 
140
            const int ya_high = qMin (mata->radius, yr + matb->radius);
 
141
            const int xa_low  = qMax (-mata->radius, xr - matb->radius);
 
142
            const int xa_high = qMin (mata->radius, xr + matb->radius);
 
143
            register double val = 0.0;
 
144
 
 
145
            for (ya = ya_low; ya <= ya_high; ++ya)
 
146
            {
 
147
                for (xa = xa_low; xa <= xa_high; ++xa)
 
148
                {
 
149
                    val += c_mat_elt (mata, xa, ya) *
 
150
                        c_mat_elt (matb, xr - xa, yr - ya);
 
151
                }
 
152
            }
 
153
 
 
154
            *c_mat_eltptr (result, xr, yr) = val;
 
155
        }
 
156
    }
 
157
}
 
158
 
 
159
void RefocusMatrix::convolve_star_mat (CMat * result, const CMat * const mata, const CMat * const matb)
 
160
{
 
161
    register int xr, yr, xa, ya;
 
162
 
 
163
    for (yr = -result->radius; yr <= result->radius; ++yr)
 
164
    {
 
165
        for (xr = -result->radius; xr <= result->radius; ++xr)
 
166
        {
 
167
            const int ya_low = qMax (-mata->radius, -matb->radius - yr);
 
168
            const int ya_high = qMin (mata->radius, matb->radius - yr);
 
169
            const int xa_low = qMax (-mata->radius, -matb->radius - xr);
 
170
            const int xa_high = qMin (mata->radius, matb->radius - xr);
 
171
            register double val = 0.0;
 
172
 
 
173
            for (ya = ya_low; ya <= ya_high; ++ya)
 
174
            {
 
175
                for (xa = xa_low; xa <= xa_high; ++xa)
 
176
                {
 
177
                    val += c_mat_elt (mata, xa, ya) *
 
178
                        c_mat_elt (matb, xr + xa, yr + ya);
 
179
                }
 
180
            }
 
181
 
 
182
            *c_mat_eltptr (result, xr, yr) = val;
 
183
        }
 
184
    }
 
185
}
 
186
 
 
187
void RefocusMatrix::convolve_mat_fun (CMat * result, const CMat * const mata, double (f) (int, int))
 
188
{
 
189
    register int xr, yr, xa, ya;
 
190
 
 
191
    for (yr = -result->radius; yr <= result->radius; ++yr)
 
192
    {
 
193
        for (xr = -result->radius; xr <= result->radius; ++xr)
 
194
        {
 
195
            register double val = 0.0;
 
196
 
 
197
            for (ya = -mata->radius; ya <= mata->radius; ++ya)
 
198
            {
 
199
                for (xa = -mata->radius; xa <= mata->radius; ++xa)
 
200
                {
 
201
                    val += c_mat_elt (mata, xa, ya) * f (xr - xa, yr - ya);
 
202
                }
 
203
            }
 
204
 
 
205
            *c_mat_eltptr (result, xr, yr) = val;
 
206
        }
 
207
    }
 
208
}
 
209
 
 
210
int RefocusMatrix::as_idx (const int k, const int l, const int m)
 
211
{
 
212
    return ((k + m) * (2 * m + 1) + (l + m));
 
213
}
 
214
 
 
215
int RefocusMatrix::as_cidx (const int k, const int l)
 
216
{
 
217
    const int a = qMax (qAbs (k), qAbs (l));
 
218
    const int b = qMin (qAbs (k), qAbs (l));
 
219
    return ((a * (a + 1)) / 2 + b);
 
220
}
 
221
 
 
222
void RefocusMatrix::print_c_mat (const CMat * const mat)
 
223
{
 
224
    register int x, y;
 
225
 
 
226
    for (y = -mat->radius; y <= mat->radius; ++y)
 
227
    {
 
228
        QString output, num;
 
229
 
 
230
        for (x = -mat->radius; x <= mat->radius; ++x)
 
231
        {
 
232
            output.append( num.setNum( c_mat_elt (mat, x, y) ) );
 
233
        }
 
234
 
 
235
        kDebug() << output;
 
236
    }
 
237
}
 
238
 
 
239
void RefocusMatrix::print_matrix (Mat * matrix)
 
240
{
 
241
    int col_idx, row_idx;
 
242
 
 
243
    for (row_idx = 0; row_idx < matrix->rows; ++row_idx)
 
244
    {
 
245
        QString output, num;
 
246
 
 
247
        for (col_idx = 0; col_idx < matrix->cols; ++col_idx)
 
248
        {
 
249
            output.append( num.setNum( mat_elt (matrix, row_idx, col_idx) ) );
 
250
        }
 
251
 
 
252
        kDebug() << output;
 
253
    }
 
254
}
 
255
 
 
256
Mat *RefocusMatrix::make_s_matrix (CMat * mat, int m, double noise_factor)
 
257
{
 
258
    const int mat_size = SQR (2 * m + 1);
 
259
    Mat *result = allocate_matrix (mat_size, mat_size);
 
260
    register int yr, yc, xr, xc;
 
261
 
 
262
    for (yr = -m; yr <= m; ++yr)
 
263
    {
 
264
        for (xr = -m; xr <= m; ++xr)
 
265
        {
 
266
            for (yc = -m; yc <= m; ++yc)
 
267
            {
 
268
                for (xc = -m; xc <= m; ++xc)
 
269
                {
 
270
                    *mat_eltptr (result, as_idx (xr, yr, m),
 
271
                                as_idx (xc, yc, m)) =
 
272
                        c_mat_elt (mat, xr - xc, yr - yc);
 
273
                    if ((xr == xc) && (yr == yc))
 
274
                    {
 
275
                        *mat_eltptr (result, as_idx (xr, yr, m),
 
276
                                    as_idx (xc, yc, m)) += noise_factor;
 
277
                    }
 
278
                }
 
279
            }
 
280
        }
 
281
    }
 
282
 
 
283
    return (result);
 
284
}
 
285
 
 
286
Mat *RefocusMatrix::make_s_cmatrix (CMat * mat, int m, double noise_factor)
 
287
{
 
288
    const int mat_size = as_cidx (m + 1, 0);
 
289
    Mat *result = allocate_matrix (mat_size, mat_size);
 
290
    register int yr, yc, xr, xc;
 
291
 
 
292
    for (yr = 0; yr <= m; ++yr)
 
293
    {
 
294
        for (xr = 0; xr <= yr; ++xr)
 
295
        {
 
296
            for (yc = -m; yc <= m; ++yc)
 
297
            {
 
298
                for (xc = -m; xc <= m; ++xc)
 
299
                {
 
300
                    *mat_eltptr (result, as_cidx (xr, yr), as_cidx (xc, yc)) +=
 
301
                        c_mat_elt (mat, xr - xc, yr - yc);
 
302
                    if ((xr == xc) && (yr == yc))
 
303
                    {
 
304
                        *mat_eltptr (result, as_cidx (xr, yr),
 
305
                                    as_cidx (xc, yc)) += noise_factor;
 
306
                    }
 
307
                }
 
308
            }
 
309
        }
 
310
    }
 
311
 
 
312
    return (result);
 
313
}
 
314
 
 
315
double RefocusMatrix::correlation (const int x, const int y, const double gamma, const double musq)
 
316
{
 
317
    return (musq + pow ((double)gamma, (double)sqrt (SQR (x) + SQR (y))));
 
318
}
 
319
 
 
320
Mat *RefocusMatrix::copy_vec (const CMat * const mat, const int m)
 
321
{
 
322
    Mat *result = allocate_matrix (SQR (2 * m + 1), 1);
 
323
    register int x, y, index = 0;
 
324
 
 
325
    for (y = -m; y <= m; ++y)
 
326
    {
 
327
        for (x = -m; x <= m; ++x)
 
328
        {
 
329
            *mat_eltptr (result, index, 0) = c_mat_elt (mat, x, y);
 
330
            ++index;
 
331
        }
 
332
    }
 
333
 
 
334
    Q_ASSERT (index == SQR (2 * m + 1));
 
335
    return (result);
 
336
}
 
337
 
 
338
Mat *RefocusMatrix::copy_cvec (const CMat * const mat, const int m)
 
339
{
 
340
    Mat *result = allocate_matrix (as_cidx (m + 1, 0), 1);
 
341
    register int x, y, index = 0;
 
342
 
 
343
    for (y = 0; y <= m; ++y)
 
344
    {
 
345
        for (x = 0; x <= y; ++x)
 
346
        {
 
347
            *mat_eltptr (result, index, 0) = c_mat_elt (mat, x, y);
 
348
            ++index;
 
349
        }
 
350
    }
 
351
 
 
352
    Q_ASSERT (index == as_cidx (m + 1, 0));
 
353
    return (result);
 
354
}
 
355
 
 
356
CMat *RefocusMatrix::copy_cvec2mat (const Mat * const cvec, const int m)
 
357
{
 
358
    CMat *result = allocate_c_mat (m);
 
359
    register int x, y;
 
360
 
 
361
    for (y = -m; y <= m; ++y)
 
362
    {
 
363
        for (x = -m; x <= m; ++x)
 
364
        {
 
365
            *c_mat_eltptr (result, x, y) = mat_elt (cvec, as_cidx (x, y), 0);
 
366
        }
 
367
    }
 
368
 
 
369
    return (result);
 
370
}
 
371
 
 
372
CMat *RefocusMatrix::copy_vec2mat (const Mat * const cvec, const int m)
 
373
{
 
374
    CMat *result = allocate_c_mat (m);
 
375
    register int x, y;
 
376
 
 
377
    for (y = -m; y <= m; ++y)
 
378
    {
 
379
        for (x = -m; x <= m; ++x)
 
380
        {
 
381
            *c_mat_eltptr (result, x, y) = mat_elt (cvec, as_idx (x, y, m), 0);
 
382
        }
 
383
    }
 
384
 
 
385
    return (result);
 
386
}
 
387
 
 
388
CMat *RefocusMatrix::compute_g (const CMat * const convolution, const int m, const double gamma,
 
389
                                const double noise_factor, const double musq, const bool symmetric)
 
390
{
 
391
    CMat h_conv_ruv, a, corr;
 
392
    CMat *result;
 
393
    Mat *b;
 
394
    Mat *s;
 
395
    int status;
 
396
 
 
397
    init_c_mat (&h_conv_ruv, 3 * m);
 
398
    fill_matrix2 (&corr, 4 * m, correlation, gamma, musq);
 
399
    convolve_mat (&h_conv_ruv, convolution, &corr);
 
400
    init_c_mat (&a, 2 * m);
 
401
    convolve_star_mat (&a, convolution, &h_conv_ruv);
 
402
 
 
403
    if (symmetric)
 
404
    {
 
405
        s = make_s_cmatrix (&a, m, noise_factor);
 
406
        b = copy_cvec (&h_conv_ruv, m);
 
407
    }
 
408
    else
 
409
    {
 
410
        s = make_s_matrix (&a, m, noise_factor);
 
411
        b = copy_vec (&h_conv_ruv, m);
 
412
    }
 
413
 
 
414
#ifdef RF_DEBUG
 
415
    kDebug() << "Convolution:";
 
416
    print_c_mat (convolution);
 
417
    kDebug() << "h_conv_ruv:";
 
418
    print_c_mat (&h_conv_ruv);
 
419
    kDebug() << "Value of s:";
 
420
    print_matrix (s);
 
421
#endif
 
422
 
 
423
    Q_ASSERT (s->cols == s->rows);
 
424
    Q_ASSERT (s->rows == b->rows);
 
425
    status = dgesv (s->rows, 1, s->data, s->rows, b->data, b->rows);
 
426
 
 
427
    if (symmetric)
 
428
    {
 
429
        result = copy_cvec2mat (b, m);
 
430
    }
 
431
    else
 
432
    {
 
433
        result = copy_vec2mat (b, m);
 
434
    }
 
435
 
 
436
#ifdef RF_DEBUG
 
437
    kDebug() << "Deconvolution:";
 
438
    print_c_mat (result);
 
439
#endif
 
440
 
 
441
    finish_c_mat (&a);
 
442
    finish_c_mat (&h_conv_ruv);
 
443
    finish_c_mat (&corr);
 
444
    finish_and_free_matrix (s);
 
445
    finish_and_free_matrix (b);
 
446
    return (result);
 
447
}
 
448
 
 
449
CMat *RefocusMatrix::compute_g_matrix (const CMat * const convolution, const int m,
 
450
                                       const double gamma, const double noise_factor,
 
451
                                       const double musq, const bool symmetric)
 
452
{
 
453
#ifdef RF_DEBUG
 
454
    kDebug() << "matrix size: " << m;
 
455
    kDebug() << "correlation: " << gamma;
 
456
    kDebug() << "noise: " << noise_factor;
 
457
#endif
 
458
 
 
459
    CMat *g = compute_g (convolution, m, gamma, noise_factor, musq, symmetric);
 
460
    int r, c;
 
461
    double sum = 0.0;
 
462
 
 
463
    /* Determine sum of array */
 
464
    for (r = -g->radius; r <= g->radius; ++r)
 
465
    {
 
466
        for (c = -g->radius; c <= g->radius; ++c)
 
467
        {
 
468
            sum += c_mat_elt (g, r, c);
 
469
        }
 
470
    }
 
471
 
 
472
    for (r = -g->radius; r <= g->radius; ++r)
 
473
    {
 
474
        for (c = -g->radius; c <= g->radius; ++c)
 
475
        {
 
476
            *c_mat_eltptr (g, r, c) /= sum;
 
477
        }
 
478
    }
 
479
 
 
480
    return (g);
 
481
}
 
482
 
 
483
void RefocusMatrix::fill_matrix (CMat * matrix, const int m,
 
484
                                 double f (const int, const int, const double),
 
485
                                 const double fun_arg)
 
486
{
 
487
    register int x, y;
 
488
    init_c_mat (matrix, m);
 
489
 
 
490
    for (y = -m; y <= m; ++y)
 
491
    {
 
492
        for (x = -m; x <= m; ++x)
 
493
        {
 
494
            *c_mat_eltptr (matrix, x, y) = f (x, y, fun_arg);
 
495
        }
 
496
    }
 
497
}
 
498
 
 
499
void RefocusMatrix::fill_matrix2 (CMat * matrix, const int m,
 
500
                                  double f (const int, const int, const double, const double),
 
501
                                  const double fun_arg1, const double fun_arg2)
 
502
{
 
503
    register int x, y;
 
504
    init_c_mat (matrix, m);
 
505
 
 
506
    for (y = -m; y <= m; ++y)
 
507
    {
 
508
        for (x = -m; x <= m; ++x)
 
509
        {
 
510
            *c_mat_eltptr (matrix, x, y) = f (x, y, fun_arg1, fun_arg2);
 
511
        }
 
512
    }
 
513
}
 
514
 
 
515
void RefocusMatrix::make_gaussian_convolution (const double gradius, CMat * convolution, const int m)
 
516
{
 
517
    register int x, y;
 
518
 
 
519
#ifdef RF_DEBUG
 
520
    kDebug() << "gauss: " << gradius;
 
521
#endif
 
522
 
 
523
    init_c_mat (convolution, m);
 
524
 
 
525
    if (SQR (gradius) <= 1.0e10F / 3.40282347e28F)
 
526
    {
 
527
        for (y = -m; y <= m; ++y)
 
528
        {
 
529
            for (x = -m; x <= m; ++x)
 
530
            {
 
531
                *c_mat_eltptr (convolution, x, y) = 0;
 
532
            }
 
533
        }
 
534
 
 
535
        *c_mat_eltptr (convolution, 0, 0) = 1;
 
536
    }
 
537
    else
 
538
    {
 
539
        const double alpha = log (2.0) / SQR (gradius);
 
540
 
 
541
        for (y = -m; y <= m; ++y)
 
542
        {
 
543
            for (x = -m; x <= m; ++x)
 
544
            {
 
545
                *c_mat_eltptr (convolution, x, y) =
 
546
                    exp (-alpha * (SQR (x) + SQR (y)));
 
547
            }
 
548
        }
 
549
    }
 
550
}
 
551
 
 
552
/** Return the integral of sqrt(radius^2 - z^2) for z = 0 to x. */
 
553
 
 
554
double RefocusMatrix::circle_integral (const double x, const double radius)
 
555
{
 
556
    if (radius == 0)
 
557
    {
 
558
        // Perhaps some epsilon must be added here.
 
559
        return (0);
 
560
    }
 
561
    else
 
562
    {
 
563
        const double sin = x / radius;
 
564
        const double sq_diff = SQR (radius) - SQR (x);
 
565
        // From a mathematical point of view the following is redundant.
 
566
        // Numerically they are not equivalent!
 
567
 
 
568
        if ((sq_diff < 0.0) || (sin < -1.0) || (sin > 1.0))
 
569
        {
 
570
            if (sin < 0)
 
571
            {
 
572
                return (-0.25 * SQR (radius) * M_PI);
 
573
            }
 
574
            else
 
575
            {
 
576
                return (0.25 * SQR (radius) * M_PI);
 
577
            }
 
578
        }
 
579
        else
 
580
        {
 
581
            return (0.5 * x * sqrt (sq_diff) + 0.5 * SQR (radius) * asin (sin));
 
582
        }
 
583
    }
 
584
}
 
585
 
 
586
double RefocusMatrix::circle_intensity (const int x, const int y, const double radius)
 
587
{
 
588
    if (radius == 0)
 
589
    {
 
590
        return (((x == 0) && (y == 0)) ? 1 : 0);
 
591
    }
 
592
    else
 
593
    {
 
594
        register double xlo = qAbs (x) - 0.5, xhi = qAbs (x) + 0.5,
 
595
            ylo = qAbs (y) - 0.5, yhi = qAbs (y) + 0.5;
 
596
        register double symmetry_factor = 1, xc1, xc2;
 
597
 
 
598
        if (xlo < 0)
 
599
        {
 
600
            xlo = 0;
 
601
            symmetry_factor *= 2;
 
602
        }
 
603
 
 
604
        if (ylo < 0)
 
605
        {
 
606
            ylo = 0;
 
607
            symmetry_factor *= 2;
 
608
        }
 
609
 
 
610
        if (SQR (xlo) + SQR (yhi) > SQR (radius))
 
611
        {
 
612
            xc1 = xlo;
 
613
        }
 
614
        else if (SQR (xhi) + SQR (yhi) > SQR (radius))
 
615
        {
 
616
            xc1 = sqrt (SQR (radius) - SQR (yhi));
 
617
        }
 
618
        else
 
619
        {
 
620
            xc1 = xhi;
 
621
        }
 
622
 
 
623
        if (SQR (xlo) + SQR (ylo) > SQR (radius))
 
624
        {
 
625
            xc2 = xlo;
 
626
        }
 
627
        else if (SQR (xhi) + SQR (ylo) > SQR (radius))
 
628
        {
 
629
            xc2 = sqrt (SQR (radius) - SQR (ylo));
 
630
        }
 
631
        else
 
632
        {
 
633
            xc2 = xhi;
 
634
        }
 
635
 
 
636
        return (((yhi - ylo) * (xc1 - xlo) +
 
637
                circle_integral (xc2, radius) - circle_integral (xc1, radius) -
 
638
                (xc2 - xc1) * ylo) * symmetry_factor / (M_PI * SQR (radius)));
 
639
    }
 
640
}
 
641
 
 
642
void RefocusMatrix::make_circle_convolution (const double radius, CMat * convolution, const int m)
 
643
{
 
644
#ifdef RF_DEBUG
 
645
    kDebug() << "radius: " << radius;
 
646
#endif
 
647
 
 
648
    fill_matrix (convolution, m, circle_intensity, radius);
 
649
}
 
650
 
 
651
int RefocusMatrix::dgesv (const int N, const int NRHS, double *A, const int lda, double *B, const int ldb)
 
652
{
 
653
    int result = 0;
 
654
    integer i_N = N, i_NHRS = NRHS, i_lda = lda, i_ldb = ldb, info;
 
655
    integer *ipiv = new integer[N];
 
656
 
 
657
    // Clapack call.
 
658
    dgesv_ (&i_N, &i_NHRS, A, &i_lda, ipiv, B, &i_ldb, &info);
 
659
 
 
660
    delete [] ipiv;
 
661
    result = info;
 
662
    return (result);
 
663
}
 
664
 
 
665
}  // namespace Digikam