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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.