"Fossies" - the Fresh Open Source Software Archive 
Member "pfstools-2.2.0/src/camera/mitsunaga99_numerical.cpp" (12 Aug 2021, 3655 Bytes) of package /linux/privat/pfstools-2.2.0.tgz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
1 /**
2 * @brief Algorithm for automatic HDR images capture.
3 *
4 *
5 * This file is a part of PFS CALIBRATION package.
6 * ----------------------------------------------------------------------
7 * Copyright (C) 2006 Radoslaw Mantiuk
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 * ----------------------------------------------------------------------
23 *
24 * @author Radoslaw Mantiuk, <radoslaw.mantiuk@gmail.com>
25 *
26 */
27
28
29 #include <math.h>
30 #include <stdlib.h>
31 #include <stdio.h>
32 #include <vector>
33 #define NRANSI
34 #include "nrutil.h"
35
36 #include "mitsunaga99_numerical.h"
37
38
39 #define TINY 1.0e-20;
40 void ludcmp(float **a, int n, int *indx, float *d)
41 {
42 int i,imax,j,k;
43 float big,dum,sum,temp;
44 float *vv;
45
46 vv=vector(1,n);
47 *d=1.0;
48 for (i=1;i<=n;i++) {
49 big=0.0;
50 for (j=1;j<=n;j++)
51 if ((temp=fabs(a[i][j])) > big) big=temp;
52 if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
53 vv[i]=1.0/big;
54 }
55
56 for (j=1;j<=n;j++) {
57 for (i=1;i<j;i++) {
58 sum=a[i][j];
59 for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
60 a[i][j]=sum;
61 }
62 big=0.0;
63 for (i=j;i<=n;i++) {
64 sum=a[i][j];
65 for (k=1;k<j;k++)
66 sum -= a[i][k]*a[k][j];
67 a[i][j]=sum;
68 if ( (dum=vv[i]*fabs(sum)) >= big) {
69 big=dum;
70 imax=i;
71 }
72 }
73 if (j != imax) {
74 for (k=1;k<=n;k++) {
75 dum=a[imax][k];
76 a[imax][k]=a[j][k];
77 a[j][k]=dum;
78 }
79 *d = -(*d);
80 vv[imax]=vv[j];
81 }
82 indx[j]=imax;
83 if (a[j][j] == 0.0) a[j][j]=TINY;
84 if (j != n) {
85 dum=1.0/(a[j][j]);
86 for (i=j+1;i<=n;i++) a[i][j] *= dum;
87 }
88 }
89 free_vector(vv,1,n);
90 }
91 #undef TINY
92
93 void lubksb(float **a, int n, int *indx, float b[])
94 {
95 int i,ii=0,ip,j;
96 float sum;
97
98 for (i=1;i<=n;i++) {
99 ip=indx[i];
100 sum=b[ip];
101 b[ip]=b[i];
102 if (ii)
103 for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
104 else if (sum) ii=i;
105 b[i]=sum;
106 }
107 for (i=n;i>=1;i--) {
108 sum=b[i];
109 for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
110 b[i]=sum/a[i][i];
111 }
112 }
113
114
115 /** Solves the system of linear equations. The result is returned in b[].
116 * @param n number of equations (equal to number of variables)
117 * @param a[][] coefficients of variables
118 * @param b[] free coefficients
119 */
120 int MitsunagaNumerical::linearEquationsSystem( int n, float a[][NR_LINEAR_EQUATIONS_MAX], float b[]) {
121
122 if( n >= NR_LINEAR_EQUATIONS_MAX) {
123 fprintf( stderr, "WARNING: too many equations (nr_solver_linear_equations())\n");
124 return 1;
125 }
126
127 float dd;
128
129 int N = n+1;
130
131 int *indx = new int[N];
132 float *bb = new float[N];
133 // int indx[N];
134 // float bb[N];
135
136 float** aa;
137 aa = (float**)malloc(sizeof(float*) * N);
138 for(int i = 0; i < N; i++)
139 aa[i] = (float*)malloc(sizeof(float) * N);
140
141 for( int i = 0; i < n; i++) {
142 bb[i+1] = b[i];
143
144 for( int j = 0; j < n; j ++) {
145 aa[i+1][j+1] = a[i][j];
146 }
147 }
148
149 ludcmp( aa, n, indx, &dd);
150 lubksb( aa, n, indx, bb);
151
152 for(int i = 0; i < N; i++)
153 if( aa[i] != NULL)
154 free(aa[i]);
155 if( aa != NULL)
156 free(aa);
157
158
159 for( int i = 0; i < n; i++) {
160 b[i] = bb[i+1];
161 }
162
163 delete[] indx;
164 delete[] bb;
165
166 return 0;
167 }
168
169
170
171