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

Linux Cross Reference
Tina4/src/vision/improc/im_poly.c

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

  1 #include <math.h>
  2 #include <tina/sys.h>
  3 #include <tina/sysfuncs.h>
  4 #include <tina/math.h>
  5 #include <tina/mathfuncs.h>
  6 #include <tina/vision.h>
  7 #include <tina/visionfuncs.h>
  8 #include <tina/hist_funcs.h>
  9 
 10 static Imrect *inim=NULL;
 11 static Imrect *outim=NULL;
 12 static shistogram *peaks=NULL;
 13 static int n_par = 4, m_par = 3;
 14 static double sigma = 1.0;
 15 
 16 shistogram *get_peaks_hist(void)
 17 {
 18     return(peaks);
 19 }
 20 
 21 double poly_xy(int n, double *a, float x, float y)
 22 {
 23     int i, count, s;
 24     double *z;
 25     double total=0.0;
 26     z = (double *)ralloc(n*sizeof(double));
 27 
 28     z[0] = 1.0;
 29     s = 1;
 30     count = 1;
 31     while (count<n)
 32     {
 33        for (i=0;i<s && count<n ;i++, count++)
 34        {
 35           z[count] = x*z[count-s];  
 36        }
 37        if (count>=n) continue;
 38        s++;
 39        z[count] = y*z[count-s];
 40        count++;
 41     }
 42     for (i=0;i<n;i++)
 43     {
 44        total += z[i]*a[i];
 45     }
 46     rfree(z);
 47     return(total);
 48 }
 49 
 50 void imf_poly(Imrect *im1, int n, double *a)
 51 {
 52     float  *row1;
 53     Imregion *roi;
 54     double x,y, cenx, ceny;
 55     int     lx, ux, ly, uy;
 56     int     i, j;
 57 
 58     roi = im1->region;
 59     lx = roi->lx;
 60     ux = roi->ux;
 61     cenx = (lx +ux )/2.0;
 62     ly = roi->ly;
 63     uy = roi->uy;
 64     ceny = (ly +uy )/2.0;
 65     row1 = fvector_alloc(lx, ux);
 66 
 67     for (i = ly; i < uy; ++i)
 68     {
 69         y = i+0.5;
 70         for (j = lx; j < ux; ++j)
 71         {
 72             x = j+0.5;
 73             row1[j] = (float)poly_xy(n, a, x-cenx, y-ceny);
 74         }
 75         im_put_rowf(row1, im1, i, lx, ux);
 76     }
 77 
 78     fvector_free((void *) row1, lx);
 79 }
 80 
 81 double trend_hfit(int n, double *a, float x)
 82 {
 83    int j;
 84    double gtemp, gtemp1, gtemp2;
 85    double dgauss(double a, double b, double c);
 86 
 87    gtemp1 = 0.0;
 88    gtemp2 = 0.0;
 89    for (j=0;j<m_par;j++)
 90    {
 91       gtemp = dgauss(a[j], sigma, x);
 92       gtemp1 += gtemp;
 93       gtemp2 += gtemp*gtemp;
 94    }
 95    if (gtemp1!=0.0) gtemp = gtemp2/gtemp1;
 96    else gtemp = 0.0;
 97    return (100*gtemp);
 98 }
 99 
100 static double pixchisq(int n, double *a)
101 {
102     double chisq = 0.0, temp;
103     Imregion *roi;
104     float  *row1, *row2;
105     int     lx, ux, ly, uy;
106     int     width, height;
107     int     i, j;
108 
109     hreset(peaks);
110     imf_poly(outim, n_par, a);
111     roi = inim->region;
112     lx = roi->lx;
113     ux = roi->ux;
114     ly = roi->ly;
115     uy = roi->uy;
116 
117     row1 = fvector_alloc(lx, ux);
118     row2 = fvector_alloc(lx, ux);
119     for (i = ly; i < uy; ++i)
120     {
121         im_get_rowf(row1, inim, i, lx, ux);
122         im_get_rowf(row2, outim, i, lx, ux);
123         for (j = lx; j < ux; ++j)
124         {
125             if (row1[j]>(-1000.0)) /* cheat to ignore pixels */
126             {
127                 hfill1s(peaks, row1[j] * row2[j], 1.0);           
128             }
129         }
130     }
131     for (i=0;i<peaks->xbins;i++)
132     {
133         temp = peaks->array[0][i];
134         chisq -= log(1.0+temp)*trend_hfit(m_par,&a[n_par],peaks->xmin+(i+0.5)*peaks->xincr);
135     }
136  
137     fvector_free((void *) row1, lx);
138     fvector_free((void *) row2, lx);
139 
140     return(chisq);
141 }
142 
143 Imrect* trend_fit(Imrect *im, int m, double *am)
144 {
145     Imrect *im1;
146     Imregion *roi=NULL;
147     double *a;
148     double c_test1 = 0.00001;
149     double *lambda;
150     int i;
151 
152     m_par = m-1;
153     sigma = am[m_par];
154     a = (double *)ralloc((n_par+m_par)*sizeof(double));
155     lambda = (double *)ralloc((n_par+m_par)*sizeof(double));
156     for (i=0;i<m_par;i++)
157     {
158        a[i+n_par] = am[i];
159        lambda[i+n_par] = 0.0;
160     }
161     lambda[n_par] = 0.0;
162 
163     if (peaks !=NULL) peaks = hfree(peaks);
164     peaks = hbook1(0,"peaks/0",a[n_par]-3.0*sigma,a[n_par+m_par-1]+3.0*sigma,128);
165     im1 = im_alloc(im->height, im->width, im->region, float_v);
166     inim = im;
167     outim = im1;
168     a[0] = 1.0;
169     a[1] = 0.0;
170     a[2] = 0.0;
171     a[3] = 0.0;
172 /*
173     a[4] = 0.0;
174     a[5] = 0.0;
175 */
176     lambda[0] = 0.0;
177     lambda[1] = 0.0001;
178     lambda[2] = 0.0001;
179     lambda[3] = 0.000001;
180 /*
181     lambda[4] = 0.000001;
182     lambda[5] = 0.000001;
183 */
184 
185     pixchisq(n_par+m_par,a);
186     simplexmin2(n_par+m_par, a, lambda, pixchisq, c_test1, (void (*) (char *)) format); 
187     pixchisq(n_par+m_par,a);
188 
189     for (i=0;i<m_par;i++)
190     {
191        am[i] = a[i+n_par];
192     }
193     hsuper(peaks, m_par, trend_hfit, am);
194     rfree(a);
195     outim=NULL;
196     return(im1);
197 }
198 

~ [ 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.