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

Linux Cross Reference
Tina4/src/tools/smartROI/sroi_search.c

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

  1 /*
  2  * sroi_search.c
  3  *
  4  * code for simplex image search
  5  *
  6  */
  7 
  8 #include <stdio.h>
  9 #include <values.h>
 10 #include <math.h>
 11 #include <tina/sys.h>
 12 #include <tina/sysfuncs.h>
 13 #include <tina/math.h>
 14 #include <tina/mathfuncs.h>
 15 #include <tina/vision.h>
 16 #include <tina/visionfuncs.h>
 17 #include <tina/tv.h>
 18 #include <tina/tvfuncs.h>
 19 #include <tina/draw.h>
 20 #include <tina/drawfuncs.h>
 21 #include <tina/infrtoolfuncs.h>
 22 #include <tina/toolsfuncs.h>
 23 #include <tina/vector.h>
 24 #include <tina/sroi_tool.h>
 25 
 26 
 27 /*
 28  * Local Variables
 29  *
 30  */
 31 
 32 static Imrect *stored_im_grad;
 33 static double *im1 = NULL, *im2 = NULL, **mean = NULL, *gim = NULL, *g_mean = NULL;
 34 static int     oldlength, oldmodelc;
 35 static int     length, ng;
 36 
 37 
 38 static void init_cfparams()
 39 {
 40     Model   *mdl;
 41 
 42     if ((mdl = get_model()) == NULL)
 43         return;
 44     if (im1 != NULL)
 45       dvector_free(im1, 0);
 46     if (im2 != NULL)
 47       dvector_free(im2, 0);
 48     if (gim != NULL)
 49       dvector_free(gim, 0);
 50     if (g_mean != NULL)
 51       dvector_free(g_mean, 0);
 52     if (mean != NULL)
 53       darray_free(mean, 0, oldmodelc, 0, oldlength);
 54 
 55     length = 2 * mdl->p_length - 1;
 56     ng = mdl->c * length;
 57     im1 = dvector_alloc(0, length);
 58     im2 = dvector_alloc(0, length);
 59     mean = darray_alloc(0, 0, mdl->c, length);
 60     gim = dvector_alloc(0, ng);
 61     g_mean = dvector_alloc(0, ng);
 62     oldlength = length;
 63     oldmodelc = mdl->c;
 64 }
 65 
 66 
 67 int model_limits()
 68 {
 69    Model *mdl = NULL;
 70 
 71    Vector **P_Eval;
 72    double  *M_Eval;
 73    double  *mw;
 74    double **pw;
 75    int      i, j;
 76 
 77    if ((mdl = get_model()) == NULL)
 78       return(0);
 79 
 80    P_Eval = mdl->P_Eval;
 81    M_Eval = (double *) mdl->M_Eval->data;
 82    mw = (double *) mdl->m_weight->data;
 83    pw = mdl->p_weight->el.double_v;
 84 
 85    if (mdl->theta < -mdl->theta_max || mdl->theta > mdl->theta_max 
 86        || mdl->s < mdl->s_min || mdl->s > mdl->s_max)
 87      return(0);
 88 
 89    for (i = 0; i < mdl->m_modes; i++)
 90    {
 91       if (mw[i] / M_Eval[i] < -mdl->m_max)
 92          return (0);
 93       /*
 94        * mw[i] = -mdl->m_max * M_Eval[i]; 
 95        */
 96       if (mw[i] / M_Eval[i] > mdl->m_max)
 97          return (0);
 98       /*
 99        * mw[i] = mdl->m_max * M_Eval[i]; 
100        */
101     }
102    return (1);
103 }
104 
105 
106 static double cost_func(int n, double *t)
107 {
108     Model   *mdl;
109     Matrix **P_Evec;
110     Vector **P_Eval;
111     double  *a, **el, **w;
112     double **g_evec, *g_eval;
113     double   cos_ai, sin_ai;
114     double   middlem1;
115     double   x0, y0, x, y, cost = 0.0, norm, **evec, *eval;
116     float  **p;
117     int      m, g, i, j, k;
118     int      start;
119 
120     if ((mdl = get_model()) == NULL)
121         return (MAXFLOAT);
122 
123     update_search_params(t, n);
124 
125     if (!model_limits())
126        return (MAXFLOAT);
127     
128     position();
129     angle_points();
130 
131     middlem1 = (double) (mdl->p_length - 1);
132     start = mdl->vsize - mdl->p_length;
133         
134     el = mdl->m->el.double_v;
135     a = (double *) mdl->alpha->data;
136     p = mdl->profile->el.float_v;
137     g_evec = mdl->G_Evec->el.double_v;
138     g_eval = (double *) mdl->G_Eval->data;
139     P_Evec = mdl->P_Evec;
140     P_Eval = mdl->P_Eval;
141     w = mdl->p_weight->el.double_v;
142 
143     if (prof_profile_get())
144     {
145         for (i = 0, g = 0; i < mdl->c; i++)
146         {
147             cos_ai = cos(a[i]);
148             sin_ai = sin(a[i]);
149             x0 = cos_ai * middlem1;
150             y0 = sin_ai * middlem1;
151             for (j = 0; j < length; j++)
152             {
153                 x = cos_ai * (double) j - x0;
154                 y = sin_ai * (double) j - y0;
155                 im1[j] = (double) im_sub_pixf(stored_im_grad,
156                                          (el[1][i] + y), (el[0][i] + x));
157 
158                 mean[i][j] = (double) p[i][start + j];
159                 im2[j] = im1[j];
160             }
161             if (prof_normalise_get())
162             {
163                 for (j = 0, norm = 0; j < length; j++)
164                     norm += fabs(im2[j]);
165                 norm /= (double) length;
166                 for (j = 0; j < length; j++)
167                     im2[j] /= norm;
168             }
169             if (prof_global_get())
170             {
171                 for (j = 0; j < length; j++)
172                 {
173                     g_mean[g] = mean[i][j];
174                     gim[g++] = im2[j];
175                 }
176             } 
177             else
178             {
179                 evec = mdl->P_Evec[i]->el.double_v;
180                 eval = (double *) mdl->P_Eval[i]->data;
181                 for (k = 0; k < mdl->p_modes; k++)
182                 {
183                     w[i][k] = 0;
184                     for (j = 0; j < length; j++)
185                     {
186                         w[i][k] += evec[j][k] * (im2[j] - mean[i][j]);
187                     }
188                     if (w[i][k] / eval[k] < -mdl->p_max)
189                         w[i][k] = -mdl->p_max * eval[k];
190                     if (w[i][k] / eval[k] > mdl->p_max)
191                         w[i][k] = mdl->p_max * eval[k];
192                 }
193                 for (j = 0; j < length; j++)
194                     for (k = 0; k < mdl->p_modes; k++)
195                         mean[i][j] += evec[j][k] * w[i][k];
196                 if (prof_mse_get())
197                     for (j = 0; j < length; j++)
198                         cost += ((im2[j] - mean[i][j]) * (im2[j] - mean[i][j]));
199                 else
200                     for (j = 0; j < length; j++)
201                         cost += fabs(im2[j] - mean[i][j]);
202             }
203         }
204         if (prof_global_get())
205         {
206             for (k = 0; k < mdl->p_modes; k++)
207             {
208                 w[0][k] = 0;
209                 for (j = 0; j < ng; j++)
210                     w[0][k] += g_evec[j][k] * (gim[j] - g_mean[j]);
211             }
212             for (k = 0; k < mdl->p_modes; k++)
213             {
214                 if (w[0][k] / g_eval[k] < -mdl->g_max)
215                     w[0][k] = -mdl->g_max * g_eval[k];
216                 if (w[0][k] / g_eval[k] > mdl->g_max)
217                     w[0][k] = mdl->g_max * g_eval[k];
218             }
219             for (j = 0; j < ng; j++)
220                 for (k = 0; k < mdl->p_modes; k++)
221                     g_mean[j] += g_evec[j][k] * w[0][k];
222             if (prof_mse_get())
223                 for (j = 0; j < ng; j++)
224                     cost += ((gim[j] - g_mean[j]) * (gim[j] - g_mean[j]));
225             else
226                 for (j = 0; j < ng; j++)
227                     cost += fabs(gim[j] - g_mean[j]);
228         }
229     } else
230     {
231         for (i = 0; i < mdl->c; i++)
232         {
233             cos_ai = cos(a[i]);
234             sin_ai = sin(a[i]);
235             x0 = -sin_ai * middlem1;
236             y0 = cos_ai * middlem1;
237             for (j = 0; j < length; j++)
238             {
239                 x = -sin_ai * (double) j - x0;
240                 y = cos_ai * (double) j - y0;
241                 cost -= (double) im_sub_pixf(stored_im_grad,
242                                          (el[1][i] + y), (el[0][i] + x));
243 
244             }
245 
246         }
247     }
248     return (cost);
249 }
250 
251 
252 int update_search_params(double *t, int size)
253 {
254     Model  *mdl = NULL;
255     double *M_Eval;
256     double *mw;
257     int     i, j;
258 
259     if ((mdl = get_model()) == NULL)
260         return(0);
261 
262     M_Eval = (double *) mdl->M_Eval->data;
263     mw = (double *) mdl->m_weight->data;
264     mdl->tx = t[0];
265     mdl->ty = t[1];
266     mdl->s = t[2];
267     mdl->theta = t[3];
268     for (i = 0; i < (size - 4) && i < mdl->m_modes; i++)
269         mw[i] = t[i + 4] * M_Eval[i];
270     for (i = (size - 4); i < mdl->m_modes; i++)
271         mw[i] = 0.0;
272 
273     return(1);
274 }
275 
276 
277 void find_simplex(Imrect *im)
278 {
279     Sroi_dparams *params = sroi_get_dparams();
280     Model   *mdl;
281     double  *x, *l, f;
282     int      n;
283 
284     x = params->sparams;
285     l = params->lambda;
286     
287     if ((mdl = get_model()) == NULL)
288       return;
289     
290     if ((n = 4 + mdl->m_modes) > PLEN)
291         return;
292 
293     if (!update_search_params(x, n))
294       return;
295         
296     if (!model_limits())
297     {
298       error("find_simplex: initial model outside limits", warning);
299       return;
300     }
301     stored_im_grad = im;
302 
303     init_cfparams();
304     f = simplexmin2(n, x, l, cost_func, params->ftol, (void (*)()) format);
305     update_search_params(x, n);
306 
307     if (!model_limits())
308       error("find_simplex: final model outside limits", warning);
309 }
310 

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