1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgPrc_poly.c,v $
37 * Date : $Date: 2004/02/13 12:30:33 $
38 * Version : $Revision: 1.7 $
39 * CVS Id : $Id: imgPrc_poly.c,v 1.7 2004/02/13 12:30:33 neil Exp $
40 *
41 * Author : Legacy TINA
42 */
43
44 /**
45 * @file
46 * @brief Fits polynomial to sparse disparity data. No longer used (NAT).
47 *
48 */
49
50 #include "imgPrc_poly.h"
51
52 #if HAVE_CONFIG_H
53 #include <config.h>
54 #endif
55
56 #include <math.h>
57 #include <tina/sys/sysDef.h>
58 #include <tina/sys/sysPro.h>
59 #include <tina/math/mathDef.h>
60 #include <tina/math/mathPro.h>
61 #include <tina/image/img_GenDef.h>
62 #include <tina/image/img_GenPro.h>
63
64 static Imrect *inim=NULL; /* static data! */
65 static Imrect *outim=NULL; /* static data! */
66 static shistogram *peaks=NULL; /* static data! */
67 static int n_par = 4, m_par = 3; /* static data! */
68 static double sigma = 1.0; /* static data! */
69
70 shistogram *get_peaks_hist(void)
71 {
72 return(peaks);
73 }
74
75 double poly_xy(int n, double *a, float x, float y)
76 {
77 int i, count, s;
78 double *z;
79 double total=0.0;
80 z = (double *)ralloc(n*sizeof(double));
81
82 z[0] = 1.0;
83 s = 1;
84 count = 1;
85 while (count<n)
86 {
87 for (i=0;i<s && count<n ;i++, count++)
88 {
89 z[count] = x*z[count-s];
90 }
91 if (count>=n) continue;
92 s++;
93 z[count] = y*z[count-s];
94 count++;
95 }
96 for (i=0;i<n;i++)
97 {
98 total += z[i]*a[i];
99 }
100 rfree(z);
101 return(total);
102 }
103
104 void imf_poly(Imrect *im1, int n, double *a)
105 {
106 float *row1;
107 Imregion *roi;
108 double x,y, cenx, ceny;
109 int lx, ux, ly, uy;
110 int i, j;
111
112 roi = im1->region;
113 lx = roi->lx;
114 ux = roi->ux;
115 cenx = (lx +ux )/2.0;
116 ly = roi->ly;
117 uy = roi->uy;
118 ceny = (ly +uy )/2.0;
119 row1 = fvector_alloc(lx, ux);
120
121 for (i = ly; i < uy; ++i)
122 {
123 y = i+0.5;
124 for (j = lx; j < ux; ++j)
125 {
126 x = j+0.5;
127 row1[j] = (float)poly_xy(n, a, x-cenx, y-ceny);
128 }
129 im_put_rowf(row1, im1, i, lx, ux);
130 }
131
132 fvector_free(row1, lx);
133 }
134
135
136 static double dgauss(double a, double b, double c)
137 {
138 double temp;
139
140 temp = (a-c)*(a-c)/(2.0*b*b);
141 return(exp(-temp));
142 }
143
144
145 double trend_hfit(int n, double *a, float x)
146 {
147 int j;
148 double gtemp, gtemp1, gtemp2;
149
150 gtemp1 = 0.0;
151 gtemp2 = 0.0;
152 for (j=0;j<m_par;j++)
153 {
154 gtemp = dgauss(a[j], sigma, x);
155 gtemp1 += gtemp;
156 gtemp2 += gtemp*gtemp;
157 }
158 if (gtemp1!=0.0) gtemp = gtemp2/gtemp1;
159 else gtemp = 0.0;
160 return (100*gtemp);
161 }
162
163 static double pixchisq(int n, double *a)
164 {
165 double chisq = 0.0, temp;
166 Imregion *roi;
167 float *row1, *row2;
168 int lx, ux, ly, uy;
169 int i, j;
170
171 hreset(peaks);
172 imf_poly(outim, n_par, a);
173 roi = inim->region;
174 lx = roi->lx;
175 ux = roi->ux;
176 ly = roi->ly;
177 uy = roi->uy;
178
179 row1 = fvector_alloc(lx, ux);
180 row2 = fvector_alloc(lx, ux);
181 for (i = ly; i < uy; ++i)
182 {
183 im_get_rowf(row1, inim, i, lx, ux);
184 im_get_rowf(row2, outim, i, lx, ux);
185 for (j = lx; j < ux; ++j)
186 {
187 if (row1[j]>(-1000.0)) /* cheat to ignore pixels */
188 {
189 hfill1s(peaks, row1[j] * row2[j], 1.0);
190 }
191 }
192 }
193 for (i=0;i<peaks->xbins;i++)
194 {
195 temp = peaks->array[0][i];
196 chisq -= log(1.0+temp)*trend_hfit(m_par,&a[n_par],peaks->xmin+(i+0.5)*peaks->xincr);
197 }
198
199 fvector_free(row1, lx);
200 fvector_free(row2, lx);
201
202 return(chisq);
203 }
204
205 Imrect* trend_fit(Imrect *im, int m, double *am)
206 {
207 Imrect *im1;
208 double *a;
209 double c_test1 = 0.00001;
210 double *lambda;
211 int i;
212
213 m_par = m-1;
214 sigma = am[m_par];
215 a = (double *)ralloc((n_par+m_par)*sizeof(double));
216 lambda = (double *)ralloc((n_par+m_par)*sizeof(double));
217 for (i=0;i<m_par;i++)
218 {
219 a[i+n_par] = am[i];
220 lambda[i+n_par] = 0.0;
221 }
222 lambda[n_par] = 0.0;
223
224 if (peaks !=NULL) peaks = hfree(peaks);
225 peaks = hbook1(0,"peaks/0",a[n_par]-3.0*sigma,a[n_par+m_par-1]+3.0*sigma,128);
226 im1 = im_alloc(im->height, im->width, im->region, float_v);
227 inim = im;
228 outim = im1;
229 a[0] = 1.0;
230 a[1] = 0.0;
231 a[2] = 0.0;
232 a[3] = 0.0;
233 /*
234 a[4] = 0.0;
235 a[5] = 0.0;
236 */
237 lambda[0] = 0.0;
238 lambda[1] = 0.0001;
239 lambda[2] = 0.0001;
240 lambda[3] = 0.000001;
241 /*
242 lambda[4] = 0.000001;
243 lambda[5] = 0.000001;
244 */
245
246 pixchisq(n_par+m_par,a);
247 simplexmin2(n_par+m_par, a, lambda, pixchisq, NULL, c_test1, (void (*) (char *)) format);
248 pixchisq(n_par+m_par,a);
249
250 for (i=0;i<m_par;i++)
251 {
252 am[i] = a[i+n_par];
253 }
254 hsuper(peaks, m_par, trend_hfit, am);
255 rfree(a);
256 outim=NULL;
257 return(im1);
258 }
259
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.