~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/image/imgPrc_poly.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.