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

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

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