1 /*******************************************************************************
2 **
3 ** Photivo
4 **
5 ** Copyright (C) 2008 Jos De Laender <jos.de_laender@telenet.be>
6 **
7 ** This file is part of Photivo.
8 **
9 ** Photivo is free software: you can redistribute it and/or modify
10 ** it under the terms of the GNU General Public License version 3
11 ** as published by the Free Software Foundation.
12 **
13 ** Photivo is distributed in the hope that it will be useful,
14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 ** GNU General Public License for more details.
17 **
18 ** You should have received a copy of the GNU General Public License
19 ** along with Photivo.  If not, see <http://www.gnu.org/licenses/>.
20 **
21 ********************************************************************************
22 **
23 ** This file is largely based on the digiKam project.
24 ** http://www.digikam.org
25 ** Copyright (C) 2005-2007 by Gilles Caulier
26 **   <caulier dot gilles at gmail dot com>
27 ** Original implementation from Refocus Gimp plug-in
28 ** Copyright (C) 1999-2003 Ernst Lippe
29 **
30 *******************************************************************************/
31 
32 #ifndef DLREFOCUSMATRIX_H
33 #define DLREFOCUSMATRIX_H
34 
35 // C ++ includes.
36 
37 #include <cstdio>
38 
39 // CMat:
40 // @radius: Radius of the matrix.
41 //
42 // Centered matrix. This is a square matrix where
43 // the indices range from [-radius, radius].
44 // The matrix contains (2 * radius + 1) ** 2 elements.
45 
46 typedef struct {
47   int     Radius;                // Radius of the matrix
48   int     RowStride;             // Size of one row = 2 * Radius + 1
49   double *Data;                  // Contents of matrix
50   double *Center;                // Points to element with index (0, 0)
51 } ptCMat;
52 
53 // Mat:
54 // @rows: Number of rows in the matrix.
55 //
56 // Normal matrix type. Indices range from
57 // [0, rows -1 ] and [0, cols - 1].
58 
59 typedef struct {
60   int     Rows;                  // Number of Rows in the matrix
61   int     Cols;                  // Number of columns in the matrix
62   double *Data;                  // Content of the matrix
63 } ptMat;
64 
65 class ptRefocusMatrix {
66 
67 public:
68 
69 static void FillMatrix(      ptCMat* Matrix,
70                        const int     MatrixRadius,
71                              double  f(int,int,double),
72                        const double  fun_arg);
73 
74 static void FillMatrix2(      ptCMat* Matrix,
75                         const int     MatrixRadius,
76                               double  f(const int,
77                                         const int,
78                                         const double,
79                                         const double),
80                         const double  fun_arg1,
81                         const double  fun_arg2);
82 
83 static void MakeCircleConvolution(const double  Radius,
84                                         ptCMat* Convolution,
85                                   const int     MatrixRadius);
86 
87 static void MakeGaussianConvolution(const double  Alpha,
88                                           ptCMat* Convolution,
89                                     const int     MatrixRadius);
90 
91 static void ConvolveStarMatrix(      ptCMat* Result,
92                                const ptCMat* const MatrixA,
93                                const ptCMat* const MatrixB);
94 
95 static ptCMat* ComputeGMatrix(const ptCMat* const Convolution,
96                               const int     MatrixRadius,
97                               const double  Gamma,
98                               const double  NoiseFactor,
99                               const double  musq,
100                               const short   Symmetric);
101 
102 static void FinishMatrix(ptMat* Matrix);
103 static void FinishAndFreeMatrix(ptMat* Matrix);
104 static void InitCMatrix(ptCMat* Matrix, const int MatrixRadius);
105 static void FinishCMatrix(ptCMat* Matrix);
106 
107 private:
108 
109 static ptMat*  AllocateMatrix(int nRows, int nCols);
110 static double* mat_eltptr(ptMat * mat, const int r, const int c);
111 static double  mat_elt(const ptMat * mat, const int r, const int c);
112 static ptCMat* AllocateCMatrix(const int Radius);
113 static inline double *c_mat_eltptr(ptCMat* mat, const int col, const int row);
114 static inline double c_mat_elt(const ptCMat* const mat,
115                                const int col,
116                                const int row);
117 static void ConvolveMatrix(      ptCMat* Result,
118                            const ptCMat* const MatrixA,
119                            const ptCMat* const MatrixB);
120 static void ConvolveMatrixFunction(      ptCMat* Result,
121                                    const ptCMat* const MatrixA,
122                                          double(f)(int, int));
123 static int as_idx(const int k, const int l, const int m);
124 static int as_cidx(const int k, const int l);
125 static ptMat *MakeSMatrix(ptCMat* Matrix, int m, double NoiseFactor);
126 static ptMat *MakeSCMatrix(ptCMat* Matrix, int m, double NoiseFactor);
127 static double Correlation(const int    x,
128                           const int    y,
129                           const double Gamma,
130                           const double musq);
131 static ptMat*  CopyVector(const ptCMat* const Matrix, const int m);
132 static ptMat*  CopyCVector(const ptCMat* const Matrix, const int m);
133 static ptCMat* CopyCVector2Matrix(const ptMat* const cvec, const int m);
134 static ptCMat* CopyVector2Matrix(const ptMat* const cvec, const int m);
135 static ptCMat* ComputeG(const ptCMat* const Convolution,
136                         const int     m,
137                         const double  Gamma,
138                         const double  NoiseFactor,
139                         const double  musq,
140                         const short   Symmetric);
141 static double CircleIntegral(const double x, const double Radius);
142 static double CircleIntensity(const int x, const int y, const double Radius);
143 // CLapack interface.
144 static int dgesv (const int     N,
145                   const int     NRHS,
146                         double* A,
147                   const int     lda,
148                         double* B,
149                   const int     ldb);
150 
151 };
152 
153 #endif
154 
155 ////////////////////////////////////////////////////////////////////////////////
156