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

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

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