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

Linux Cross Reference
Tina6/tina-libs/tina/medical/medSroi_search.c

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/medical/medSroi_search.c,v $
 27  * Date    :  $Date: 2004/03/10 14:56:22 $
 28  * Version :  $Revision: 1.7 $
 29  * CVS Id  :  $Id: medSroi_search.c,v 1.7 2004/03/10 14:56:22 neil Exp $
 30  *
 31  * Author  : Legacy TINA converted by NAT
 32  *
 33 */
 34 /**
 35 *   @file 
 36 *   @brief Routines for the location of an active shape model.
 37 *          The location process supports several basic options. Selection of the profile_on
 38 *          option will switch between a snake like potantial optimisation only at sample points 
 39 *          (use an appropriate search image) to grey level minimisation along the profile
 40 *          on the basis of an active shape model with a specified number of control parameters. 
 41 *          Use of the mse_on option selects a robust greylevel difference measure rather than
 42 *          least squares. Methods are currently not robust enough for use without subsequent
 43 *          visual confirmation of the location. The methods are therefore intended as a
 44 *          supplement to user interfaces and teaching.
 45 */
 46 
 47 #if HAVE_CONFIG_H
 48 #   include <config.h>
 49 #endif
 50 
 51 #include "medSroi_search.h"
 52 
 53 #include <float.h>
 54 #include <math.h>
 55 #include <tina/sys/sysDef.h>
 56 #include <tina/sys/sysPro.h>
 57 #include <tina/math/mathDef.h>
 58 #include <tina/math/mathPro.h>
 59 #include <tina/image/imgDef.h>
 60 #include <tina/image/imgPro.h>
 61 
 62 
 63 #include "med_SroiDef.h"
 64 #include "medSroi_pca.h"
 65 #include "medSroi_sample.h"
 66 
 67 void prof_profile_on(Sroi_dparams *params, Bool on)
 68 {
 69     params->profile_on = on;
 70 }
 71 
 72 Bool prof_profile_get(Sroi_dparams *params)
 73 {
 74     return (params->profile_on);
 75 }
 76 
 77 void prof_mse_on(Sroi_dparams *params, Bool on)
 78 {
 79     params->mse_on = on;
 80 }
 81 
 82 Bool prof_mse_get(Sroi_dparams *params)
 83 {
 84     return (params->mse_on);
 85 }
 86 
 87 int model_limits(Model *mdl)
 88 {
 89    Vector **P_Eval;
 90    double  *M_Eval;
 91    double  *mw;
 92    double **pw;
 93    int      i;
 94 
 95    if (mdl == NULL)
 96       return(0);
 97 
 98    P_Eval = mdl->P_Eval;
 99    M_Eval = (double *) mdl->M_Eval->data;
100    mw = (double *) mdl->m_weight->data;
101    pw = mdl->p_weight->el.double_v;
102 
103    if (mdl->theta < mdl->theta_min || mdl->theta > mdl->theta_max 
104        || mdl->s < mdl->s_min || mdl->s > mdl->s_max)
105      return(0);
106 
107    for (i = 0; i < mdl->m_modes; i++)
108    {
109       if (mw[i] / M_Eval[i] < -mdl->m_max)
110          return (0);
111       /*
112        * mw[i] = -mdl->m_max * M_Eval[i]; 
113        */
114       if (mw[i] / M_Eval[i] > mdl->m_max)
115          return (0);
116       /*
117        * mw[i] = mdl->m_max * M_Eval[i]; 
118        */
119     }
120    return (1);
121 }
122 
123 
124 static double cost_func(int n, double *t, void *data)
125 {
126     Model   *mdl;
127     Sroi_dparams *params;
128     Matrix **P_Evec;
129     Vector **P_Eval;
130     double  *a, **el, **w;
131     double **g_evec, *g_eval;
132     double   cos_ai, sin_ai;
133     double   middlem1;
134     double   x0, y0, x, y, cost = 0.0, norm, **evec, *eval;
135     double *im1 = NULL, *im2 = NULL, **mean = NULL, *gim = NULL, *g_mean = NULL;
136     static int     length, ng;
137     float  **p;
138     int      g, i, j, k;
139     int      start;
140 
141     if ((mdl = (Model *)data) == NULL)
142         return (FLT_MAX);
143     params = mdl->params;
144 
145     update_search_params(t, n, mdl);
146 
147     if (!model_limits(mdl))
148        return (FLT_MAX);
149     
150     asm_position(mdl);
151     asm_angle_points(mdl);
152 
153     length = 2 * mdl->pcavsize - 1;
154     ng = mdl->c * length;
155     im1 = dvector_alloc(0, length);
156     im2 = dvector_alloc(0, length);
157     mean = darray_alloc(0, 0, mdl->c, length);
158     if (prof_global_get(mdl))
159     {
160        gim = dvector_alloc(0, ng);
161        g_mean = dvector_alloc(0, ng);
162     }
163 
164 /*
165     middlem1 = (double) (mdl->p_length - 1);
166     start = mdl->vsize - mdl->p_length;
167 */
168     middlem1 = (double) (mdl->pcavsize - 1);
169     start = mdl->vsize - mdl->pcavsize;
170         
171     el = mdl->m->el.double_v;
172     a = (double *) mdl->alpha->data;
173     p = mdl->profile->el.float_v;
174     g_evec = mdl->G_Evec->el.double_v;
175     g_eval = (double *) mdl->G_Eval->data;
176     P_Evec = mdl->P_Evec;
177     P_Eval = mdl->P_Eval;
178     w = mdl->p_weight->el.double_v;
179 
180     if (prof_profile_get(params))
181     {
182         for (i = 0, g = 0; i < mdl->c; i++)
183         {
184             cos_ai = cos(a[i]);
185             sin_ai = sin(a[i]);
186             x0 = cos_ai * middlem1;
187             y0 = sin_ai * middlem1;
188             for (j = 0; j < length; j++)
189             {
190                 x = cos_ai * (double) j - x0;
191                 y = sin_ai * (double) j - y0;
192                 im1[j] = (double) im_sub_pixf(mdl->im_grad,
193                                          (el[1][i] + y), (el[0][i] + x));
194 
195                 mean[i][j] = (double) p[i][start + j];
196                 im2[j] = im1[j];
197             }
198             if (prof_normalise_get(mdl))
199             {
200                 for (j = 0, norm = 0; j < length; j++)
201                     norm += fabs(im2[j]);
202                 norm /= (double) length;
203                 for (j = 0; j < length; j++)
204                     im2[j] /= norm;
205             }
206             if (prof_global_get(mdl))
207             {
208                 for (j = 0; j < length; j++)
209                 {
210                     g_mean[g] = mean[i][j];
211                     gim[g++] = im2[j];
212                 }
213             } 
214             else
215             {
216                 evec = mdl->P_Evec[i]->el.double_v;
217                 eval = (double *) mdl->P_Eval[i]->data;
218                 for (k = 0; k < mdl->p_modes; k++)
219                 {
220                     w[i][k] = 0;
221                     for (j = 0; j < length; j++)
222                     {
223                         w[i][k] += evec[j][k] * (im2[j] - mean[i][j]);
224                     }
225 /* now test for limit of model and set to limit if beyond */
226                     if (w[i][k] / eval[k] < -mdl->p_max)
227                         w[i][k] = -mdl->p_max * eval[k];
228                     if (w[i][k] / eval[k] > mdl->p_max)
229                         w[i][k] = mdl->p_max * eval[k];
230                 }
231                 for (j = 0; j < length; j++)
232                     for (k = 0; k < mdl->p_modes; k++)
233                         mean[i][j] += evec[j][k] * w[i][k];
234                 if (prof_mse_get(params))
235                     for (j = 0; j < length; j++)
236                         cost += ((im2[j] - mean[i][j]) * (im2[j] - mean[i][j]));
237                 else
238                     for (j = 0; j < length; j++)
239                         cost += fabs(im2[j] - mean[i][j]);
240             }
241         }
242         if (prof_global_get(mdl))
243         {
244             for (k = 0; k < mdl->p_modes; k++)
245             {
246                 w[0][k] = 0;
247                 for (j = 0; j < ng; j++)
248                     w[0][k] += g_evec[j][k] * (gim[j] - g_mean[j]);
249             }
250             for (k = 0; k < mdl->p_modes; k++)
251             {
252                 if (w[0][k] / g_eval[k] < -mdl->g_max)
253                     w[0][k] = -mdl->g_max * g_eval[k];
254                 if (w[0][k] / g_eval[k] > mdl->g_max)
255                     w[0][k] = mdl->g_max * g_eval[k];
256             }
257             for (j = 0; j < ng; j++)
258                 for (k = 0; k < mdl->p_modes; k++)
259                     g_mean[j] += g_evec[j][k] * w[0][k];
260             if (prof_mse_get(params))
261                 for (j = 0; j < ng; j++)
262                     cost += ((gim[j] - g_mean[j]) * (gim[j] - g_mean[j]));
263             else
264                 for (j = 0; j < ng; j++)
265                     cost += fabs(gim[j] - g_mean[j]);
266         }
267     } else
268     {
269         for (i = 0; i < mdl->c; i++)
270         {
271             cos_ai = cos(a[i]);
272             sin_ai = sin(a[i]);
273             x0 = -sin_ai * middlem1;
274             y0 = cos_ai * middlem1;
275             for (j = 0; j < length; j++)
276             {
277                 x = -sin_ai * (double) j - x0;
278                 y = cos_ai * (double) j - y0;
279                 cost -= (double) im_sub_pixf(mdl->im_grad,
280                                          (el[1][i] + y), (el[0][i] + x));
281 
282             }
283 
284         }
285     }
286     dvector_free(im1, 0);
287     dvector_free(im2, 0);
288     darray_free(mean, 0, mdl->c, 0, length);
289     if (mean != NULL)
290     {
291         dvector_free(g_mean, 0);
292         dvector_free(gim, 0);
293     }
294     return (cost);
295 }
296 
297 
298 int update_search_params(double *t, int size, Model *mdl)
299 {
300     double *M_Eval;
301     double *mw;
302     int     i;
303 
304     if (mdl == NULL)
305         return(0);
306 
307     M_Eval = (double *) mdl->M_Eval->data;
308     mw = (double *) mdl->m_weight->data;
309     mdl->tx = t[0];
310     mdl->ty = t[1];
311     mdl->s = t[2];
312     mdl->theta = t[3];
313     for (i = 0; i < (size - 4) && i < mdl->m_modes; i++)
314         mw[i] = t[i + 4] * M_Eval[i];
315     for (i = (size - 4); i < mdl->m_modes; i++)
316         mw[i] = 0.0;
317 
318     return(1);
319 }
320 
321 
322 void find_simplex(Imrect *im, Sroi_dparams *params, Model *mdl)
323 {
324     double  *x, *l, f;
325     int      n;
326 
327     x = params->sparams;
328     l = params->lambda;
329     
330     if (mdl == NULL)
331       return;
332     mdl->params = params;
333     
334     if ((n = 4 + mdl->m_modes) > PLEN)
335         return;
336 
337     if (!update_search_params(x, n, mdl))
338       return;
339         
340     if (!model_limits(mdl))
341     {
342       error("find_simplex: initial model outside limits", warning);
343       return;
344     }
345     mdl->im_grad = im;
346 
347     f = simplexmin2(n, x, l, cost_func, mdl, params->ftol, (void (*)()) format);
348     update_search_params(x, n, mdl);
349     mdl->params = NULL;
350     mdl->im_grad = NULL;
351 
352     if (!model_limits(mdl))
353       error("find_simplex: final model outside limits", warning);
354 }
355 

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