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_pca.c,v $
23 * Date : $Date: 2006/04/11 14:15:30 $
24 * Version : $Revision: 1.10 $
25 * CVS Id : $Id: medSroi_pca.c,v 1.10 2006/04/11 14:15:30 neil Exp $
26 *
27 * Author : Legacy TINA converted by NAT
28 *
29 */
30 /**
31 * @file
32 * @brief Routines for the generation of simple active shape models intended for medical use.
33 * Supports data for separate and global shape and grey level profiles taken at
34 * fixed ample locations around inner and outer boundaries of a simple shape which
35 * is rotated and scaled to a baseline.
36 *
37 *********
38 */
39
40 #if HAVE_CONFIG_H
41 # include <config.h>
42 #endif
43
44 #include "medSroi_pca.h"
45
46 #include <math.h>
47 #include <string.h>
48 #include <tina/medical/med_SroiDef.h>
49 #include <tina/medical/medSroi_alloc.h>
50 #include <tina/medical/medSroi_io.h>
51
52 #include <tina/math/mathPro.h>
53 #include <tina/image/imgDef.h>
54 #include <tina/image/imgPro.h>
55
56 void prof_normalise_on(Model *mdl, Bool on)
57 {
58 mdl->normalise_on = on;
59 }
60
61 void prof_global_on(Model *mdl, Bool on)
62 {
63 mdl->global_on = on;
64 }
65
66 Bool prof_normalise_get(Model *mdl)
67 {
68 return (mdl->normalise_on);
69 }
70
71 Bool prof_global_get(Model *mdl)
72 {
73 return (mdl->global_on);
74 }
75
76 static void matrix_swap_all_cols(Matrix * mat)
77 {
78 int n = mat->n;
79 int i, j;
80
81 for (i = 0, j = n - 1; i < j; i++, j--)
82 matrix_swap_cols(mat, i, j);
83 }
84
85
86 static void vector_swap(Vector * vec)
87 {
88 int i, j, n = vec->n;
89 double temp;
90
91 for (i = 0, j = n - 1; i < j; i++, j--)
92 {
93 temp = vector_getf(vec, i);
94 vector_putf(vector_getf(vec, j), vec, i);
95 vector_putf(temp, vec, j);
96 }
97 }
98
99
100 static void dvector_sqrt(Vector * vec)
101 {
102 int i, n = vec->n;
103 double *el = (double *) vec->data;
104
105 for (i = 0; i < n; i++)
106 el[i] = sqrt(fabs(el[i]));
107 }
108
109
110 void input_pca(char *pcafn, Model *model)
111 {
112 FILE *f;
113 Bool global_on, normalise_on;
114 char temps[256];
115 char dtype[256];
116 int i;
117
118 if (model == NULL)
119 {
120 error("input_pca: no model present", non_fatal);
121 return;
122 }
123 if ((f = fopen(pcafn, "r")) == NULL)
124 {
125 string_append(temps, "can't open file ", pcafn, 0);
126 error(temps, non_fatal);
127 return;
128 }
129 dmatrix_scanf(f, model->M_Evec);
130 dvector_scanf(f, model->M_Eval);
131 vector_format(model->M_Eval);
132
133 fscanf(f, "%s %d", dtype, &(normalise_on));
134 if (strcmp(dtype, "normalise") != 0)
135 {
136 error("Wrong data type in file", non_fatal);
137 return;
138 }
139 model->normalise_on = normalise_on;
140
141 fscanf(f, "%s %d", dtype, &(global_on));
142 if (strcmp(dtype, "global") != 0)
143 {
144 error("Wrong data type in file", non_fatal);
145 return;
146 }
147 model->global_on = global_on;
148
149 if (global_on)
150 {
151 dmatrix_scanf(f, model->G_Evec);
152 dvector_scanf(f, model->G_Eval);
153 vector_format(model->G_Eval);
154 } else
155 {
156 for (i = 0; i < model->c; i++)
157 {
158 dmatrix_scanf(f, model->P_Evec[i]);
159 dvector_scanf(f, model->P_Eval[i]);
160 }
161 }
162 fclose(f);
163 }
164
165 void update_model_sparams(Model *mdl, Sroi_dparams *params)
166 {
167 mdl->pcavsize = params->pcavsize;
168 mdl->m_modes = params->m_modes;
169 mdl->p_modes = params->p_modes;
170 mdl->p_length =params->p_length;
171 mdl->s_min = params->s_min;
172 mdl->s_max = params->s_max;
173 mdl->theta_min = params->theta_min;
174 mdl->theta_max = params->theta_max;
175 mdl->m_max = params->m_max;
176 mdl->p_max = params->p_max;
177 mdl->g_max = params->g_max;
178 }
179
180 /*
181 * original profile length is: length
182 * used profile length is : plength
183 * P_mean is of length, P_cov is of plength
184 */
185 /**
186 * @brief Function to construct an active shapoe model from a set of sample file and generate
187 * the output model as files.
188 * @param list file is modelfn, pca parameters are stored in pcafn, mean model is stored in meanfn.
189 *
190 * The model building process supports two basic options; the choice to build a model with
191 * separated (independant) shape and grey level behaviour or a global model.
192 * or the choice to normalise the grey level profile before starting.
193 * The resulting eigen vector model must be considered as a hyperplane fit on the assumtion of
194 * independant and homogenous errors on all data values.
195 * These two modes are therefore incompatable as the global model demands that the grey level data
196 * has already been scaled so that the statistical accuracy of position and grey level are matched.
197 *
198 */
199 void pca_model(char *modelfn, char *pcafn, char *meanfn, Sroi_dparams *params)
200 {
201 Model *model=NULL;
202 Matrix *M_cov=NULL, *M_Evec=NULL, **P_cov=NULL, *P_Evec=NULL, *G_cov=NULL, *G_Evec=NULL, *R_Evec=NULL;
203 Vector *M_mean=NULL, *M_Eval=NULL, *M_xy=NULL, **P_mean=NULL, *P_Eval=NULL, *G_mean=NULL, *G_Eval=NULL,
204 *R_Eval=NULL;
205 FILE *f=NULL;
206 double tx = 0, ty = 0, s = 0, theta = 0, n = 0.0;
207 double **m=NULL, **mc=NULL, *mm=NULL, *mxy=NULL, ***pc=NULL, **pm=NULL, **gc=NULL, *gm=NULL, **pg=NULL, norm;
208 float **p=NULL;
209 char *temps=NULL, fname[512];
210 int i = 0, j, r, c, length=0, plength=0, start=0, ng=0, NX=0, NXY=0;
211 int blres;
212
213 if ((modelfn == NULL) || (pcafn == NULL) || (meanfn == NULL))
214 {
215 format("pca_model: insufficient file names provided \n");
216 return;
217 }
218
219 sprintf(fname, "%s%s", modelfn, ".blt");
220 if ((n = (double)sroi_open_buildlist(fname)) == 0)
221 {
222 format("pca_model: no sample files found \n");
223 return;
224 }
225
226 while((blres = sroi_getnext_build(fname)) > 0)
227 {
228 model = input_model(fname, NULL);
229 if (model == NULL)
230 return;
231 update_model_sparams(model, params);
232
233 m = model->m->el.double_v;
234 p = model->profile->el.float_v;
235 pg = model->g_profile->el.double_v;
236
237 /* The indented lines below were enclosed in a test i==0, but i is set to zero
238 on initialisation, is not static, and is not reset before these lines. I have
239 therefore removed the test to avoid warnings that the variables below could be
240 used uninitialised PAB */
241 if (i==0)
242 {
243 length = 2 * model->vsize - 1;
244 plength = 2 * model->pcavsize - 1;
245 start = model->vsize - model->pcavsize;
246 NX = model->c;
247 NXY = 2 * NX;
248 ng = NX * length;
249 M_mean = vector_alloc(NXY, double_v);
250 M_cov = matrix_alloc(NXY, NXY, matrix_full, double_v);
251 M_Evec = matrix_alloc(NXY, NXY, matrix_full, double_v);
252 M_Eval = vector_alloc(NXY, double_v);
253 M_xy = vector_alloc(NXY, double_v);
254 mm = (double *) M_mean->data;
255 mc = M_cov->el.double_v;
256 mxy = (double *) M_xy->data;
257
258 G_mean = vector_alloc(ng, double_v);
259 G_cov = matrix_alloc(ng, ng, matrix_full, double_v);
260 G_Evec = matrix_alloc(ng, ng, matrix_full, double_v);
261 G_Eval = vector_alloc(ng, double_v);
262 R_Evec = matrix_alloc(ng, NX, matrix_full, double_v);
263 R_Eval = vector_alloc(NX, double_v);
264 gm = (double *) G_mean->data;
265 gc = G_cov->el.double_v;
266
267 P_mean = (Vector **) pvector_alloc(0, model->c);
268 P_cov = (Matrix **) pvector_alloc(0, model->c);
269 pm = (double **) pvector_alloc(0, model->c);
270 pc = (double ***) pvector_alloc(0, model->c);
271 for (r = 0; r < model->c; r++)
272 {
273 P_mean[r] = vector_alloc(length, double_v);
274 pm[r] = (double *) P_mean[r]->data;
275 P_cov[r] = matrix_alloc(plength, plength,
276 matrix_full, double_v);
277 pc[r] = P_cov[r]->el.double_v;
278 }
279 P_Evec = matrix_alloc(plength, plength,
280 matrix_full, double_v);
281 P_Eval = vector_alloc(plength, double_v);
282 }
283
284 for (r = 0; r < NX; r++)
285 {
286 mxy[r] = m[3][r];
287 mxy[NX + r] = m[2][r];
288 }
289 for (r = 0; r < NXY; r++)
290 mm[r] += mxy[r];
291 for (r = 0; r < NXY; r++)
292 {
293 mc[r][r] += mxy[r] * mxy[r];
294 for (c = r + 1; c < NXY; c++)
295 mc[r][c] += mxy[r] * mxy[c];
296 }
297
298 if (params->gradient_on)
299 {
300 for (r = 0; r < model->c; r++)
301 {
302 pg[r][0] = (double) (p[r][1] - p[r][0]);
303 pg[r][1] = (double) (p[r][2] - p[r][0]);
304 for (c = 2; c < length - 2; c++)
305 pg[r][c] = (double) (p[r][c + 2] - p[r][c - 2]);
306 pg[r][length - 2] = (double) (p[r][length - 1] - p[r][length - 3]);
307 pg[r][length - 1] = (double) (p[r][length - 1] - p[r][length - 2]);
308 }
309 } else
310 {
311 for (r = 0; r < model->c; r++)
312 for (c = 0; c < length; c++)
313 pg[r][c] = (double) p[r][c];
314 }
315
316 if (model->normalise_on)
317 {
318 for (j = 0; j < model->c; j++)
319 {
320 for (r = 0, norm = 0; r < length; r++)
321 norm += fabs(pg[j][r]);
322 norm /= (double) length;
323 for (r = 0; r < length; r++)
324 pg[j][r] /= norm;
325 }
326 }
327
328 if (model->global_on)
329 {
330 for (r = 0, j = 0; r < model->c; r++)
331 for (c = 0; c < length; c++)
332 gm[j++] += pg[r][c];
333 for (r = 0; r < ng; r++)
334 {
335 gc[r][r] += gm[r] * gm[r];
336 for (c = r + 1; c < ng; c++)
337 gc[r][c] += gm[r] * gm[c];
338 }
339 } else
340 {
341 for (j = 0; j < model->c; j++)
342 for (r = 0; r < length; r++)
343 pm[j][r] += pg[j][r];
344 for (j = 0; j < model->c; j++)
345 {
346 for (r = 0; r < plength; r++)
347 {
348 pc[j][r][r] += pg[j][start + r] * pg[j][start + r];
349 for (c = r + 1; c < plength; c++)
350 pc[j][r][c] += pg[j][start + r] * pg[j][start + c];
351 }
352 }
353 }
354 tx += model->tx;
355 ty += model->ty;
356 s += model->s;
357 theta += model->theta;
358 printf("%d %s\n", i++, fname);
359 }
360
361 sroi_close_buildlist();
362
363 if (blres != 0)
364 {
365 error("pca_model: buildlist file not finished - aborting", warning);
366 return;
367 }
368
369 for (r = 0; r < NXY; r++)
370 mm[r] /= n;
371 for (r = 0; r < NXY; r++)
372 {
373 for (c = r; c < NXY; c++)
374 mc[r][c] = mc[r][c] / n - mm[r] * mm[c];
375 for (c = r + 1; c < NXY; c++)
376 mc[c][r] = mc[r][c];
377 }
378 for (r = 0; r < NX; r++)
379 {
380 m[3][r] = mm[r];
381 m[2][r] = mm[NX + r];
382 }
383
384 if (model->global_on)
385 {
386 for (r = 0; r < ng; r++)
387 gm[r] /= n;
388 for (r = 0; r < ng; r++)
389 {
390 for (c = r; c < ng; c++)
391 gc[r][c] = gc[r][c] / n - gm[r] * gm[c];
392 for (c = r + 1; c < ng; c++)
393 gc[c][r] = gc[r][c];
394 }
395 for (r = 0, j = 0; r < NX; r++)
396 {
397 for (c = 0; c < length; c++)
398 p[r][c] = (int) gm[j++];
399 }
400 } else
401 {
402 for (j = 0; j < model->c; j++)
403 for (r = 0; r < length; r++)
404 pm[j][r] /= n;
405 for (j = 0; j < model->c; j++)
406 {
407 for (r = 0; r < plength; r++)
408 {
409 for (c = r; c < plength; c++)
410 pc[j][r][c] = pc[j][r][c] / n - pm[j][start + r] * pm[j][start + c];
411 for (c = r + 1; c < plength; c++)
412 pc[j][c][r] = pc[j][r][c];
413 }
414 }
415 for (r = 0; r < model->c; r++)
416 {
417 for (c = 0; c < length; c++)
418 p[r][c] = (float) pm[r][c];
419 }
420 }
421
422 model->tx = tx / n;
423 model->ty = ty / n;
424 model->s = s / n;
425 model->theta = theta / n;
426
427 output_model(meanfn,model);
428
429 if ((f = fopen(pcafn, "wt")) == NULL)
430 {
431 string_append(temps, "can't open file", pcafn, 0);
432 error(temps, non_fatal);
433 return;
434 }
435 matrix_eigen_sym(M_cov, M_Eval, M_Evec);
436
437 matrix_swap_all_cols(M_Evec);
438 vector_swap(M_Eval);
439 dvector_sqrt(M_Eval);
440 vector_format(M_Eval);
441 dmatrix_print(f, M_Evec);
442 dvector_printf(f, M_Eval);
443
444 fprintf(f, "normalise %d \n", (int)model->normalise_on);
445 fprintf(f, "global %d \n", (int)model->global_on);
446 if (model->global_on)
447 {
448 matrix_eigen_sym(G_cov, G_Eval, G_Evec);
449 matrix_swap_all_cols(G_Evec);
450 vector_swap(G_Eval);
451 dvector_sqrt(G_Eval);
452 vector_format(G_Eval);
453 for (c = 0; c < model->c; c++)
454 {
455 for (r = 0; r < ng; r++)
456 R_Evec->el.double_v[r][c] = G_Evec->el.double_v[r][c];
457 ((double *) (R_Eval->data))[c] = ((double *) (G_Eval->data))[c];
458 }
459 dmatrix_print(f, R_Evec);
460 dvector_printf(f, R_Eval);
461 } else
462 {
463 for (j = 0; j < model->c; j++)
464 {
465 matrix_eigen_sym(P_cov[j], P_Eval, P_Evec);
466 matrix_swap_all_cols(P_Evec);
467 vector_swap(P_Eval);
468 dvector_sqrt(P_Eval);
469 dmatrix_print(f, P_Evec);
470 dvector_printf(f, P_Eval);
471 }
472 }
473 vector_free(M_mean);
474 matrix_free(M_cov);
475 vector_free(M_Eval);
476 matrix_free(M_Evec);
477 vector_free(G_mean);
478 matrix_free(G_cov);
479 vector_free(G_Eval);
480 matrix_free(G_Evec);
481 vector_free(R_Eval);
482 matrix_free(R_Evec);
483 vector_free(P_Eval);
484 matrix_free(P_Evec);
485 for (i = 0; i < model->c; i++)
486 {
487 vector_free(P_mean[i]);
488 matrix_free(P_cov[i]);
489 }
490 fclose(f);
491 }
492
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.