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

Linux Cross Reference
Tina5/tina-libs/tina/medical/medSroi_pca.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_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 

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