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

Linux Cross Reference
Tina4/src/tools/smartROI/sroi_pca.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /*
  2  * sroi_pca.c
  3  *
  4  * pca code
  5  *
  6  */
  7 
  8 #include <stdio.h>
  9 #include <values.h>
 10 #include <math.h>
 11 #include <tina/sys.h>
 12 #include <tina/sysfuncs.h>
 13 #include <tina/math.h>
 14 #include <tina/mathfuncs.h>
 15 #include <tina/vision.h>
 16 #include <tina/visionfuncs.h>
 17 #include <tina/tv.h>
 18 #include <tina/tvfuncs.h>
 19 #include <tina/draw.h>
 20 #include <tina/drawfuncs.h>
 21 #include <tina/infrtoolfuncs.h>
 22 #include <tina/toolsfuncs.h>
 23 #include <tina/vector.h>
 24 #include <tina/sroi_tool.h>
 25 
 26 static Bool gradient_on;
 27 static Bool normalise_on;
 28 static Bool global_on;
 29 static Bool mse_on;
 30 static Bool profile_on;
 31 
 32 void prof_profile_on(Bool on)
 33 {
 34     profile_on = on;
 35 }
 36 
 37 void prof_gradient_on(Bool on)
 38 {
 39     gradient_on = on;
 40 }
 41 
 42 void prof_normalise_on(Bool on)
 43 {
 44     normalise_on = on;
 45 }
 46 
 47 void prof_mse_on(Bool on)
 48 {
 49     mse_on = on;
 50 }
 51 
 52 void prof_global_on(Bool on)
 53 {
 54     global_on = on;
 55 }
 56 
 57 Bool prof_profile_get()
 58 {
 59     return (profile_on);
 60 }
 61 
 62 Bool prof_gradient_get()
 63 {
 64     return (gradient_on);
 65 }
 66 
 67 Bool prof_normalise_get()
 68 {
 69     return (normalise_on);
 70 }
 71 
 72 Bool prof_mse_get()
 73 {
 74     return (mse_on);
 75 }
 76 
 77 Bool prof_global_get()
 78 {
 79     return (global_on);
 80 }
 81 
 82 Matrix *matrix_swap_all_cols(Matrix * mat)
 83 {
 84     int n = mat->n;
 85     int i, j;
 86 
 87     for (i = 0, j = n - 1; i < j; i++, j--)
 88         matrix_swap_cols(mat, i, j);
 89 }
 90 
 91 
 92 Vector *vector_swap(Vector * vec)
 93 {
 94     int i, j, n = vec->n;
 95     double temp;
 96 
 97     for (i = 0, j = n - 1; i < j; i++, j--)
 98     {
 99         temp = vector_getf(vec, i);
100         vector_putf(vector_getf(vec, j), vec, i);
101         vector_putf(temp, vec, j);
102     }
103 }
104 
105 
106 void dvector_sqrt(Vector * vec)
107 {
108     int i, n = vec->n;
109     double *el = (double *) vec->data;
110 
111     for (i = 0; i < n; i++)
112         el[i] = sqrt(fabs(el[i]));
113 }
114 
115 
116 void input_pca(char *pcafn)
117 {
118     Model *model;
119     FILE  *f;
120     char   temps[256];
121     int    i;
122 
123     if ((model = get_model()) == NULL)
124     {
125         error("input_pca: no model present", non_fatal);
126         return;
127     }
128     if ((f = fopen(pcafn, "r")) == NULL)
129     {
130         string_append(temps, "can't open file ", pcafn, 0);
131         error(temps, non_fatal);
132         return;
133     }
134     dmatrix_scanf(f, model->M_Evec);
135     dvector_scanf(f, model->M_Eval);
136     vector_format(model->M_Eval);
137     if (global_on)
138     {
139         dmatrix_scanf(f, model->G_Evec);
140         dvector_scanf(f, model->G_Eval);
141         vector_format(model->G_Eval);
142     } else
143     {
144         for (i = 0; i < model->c; i++)
145         {
146             dmatrix_scanf(f, model->P_Evec[i]);
147             dvector_scanf(f, model->P_Eval[i]);
148         }
149     }
150     fclose(f);
151 }
152 
153 
154 /*
155  * original profile length is: length
156  * used profile length is    : plength
157  * P_mean is of length, P_cov is of plength
158  */
159 void pca_model(char *modelfn, char *pcafn, char *meanfn)
160 {
161     Model    *model;
162     Matrix   *M_cov, *M_Evec, **P_cov, *P_Evec, *G_cov, *G_Evec, *R_Evec;
163     Vector   *M_mean, *M_Eval, *M_xy, **P_mean, *P_Eval, *G_mean, *G_Eval,
164              *R_Eval;
165     FILE     *f;
166     double    tx = 0, ty = 0, s = 0, theta = 0, n = 0.0;
167     double  **m, **mc, *mm, *mxy, ***pc, **pm, **gc, *gm, **pg, norm;
168     float   **p;
169     char     *temps, fname[512];
170     int       i = 0, j, r, c, length, plength, start, ng, NX, NXY;
171     int       blres;
172 
173     if ((modelfn == NULL) || (pcafn == NULL) || (meanfn == NULL))
174       return;
175 
176     sprintf(fname, "%s%s", modelfn, ".blt");
177     if ((n = (double)sroi_open_buildlist(fname)) == 0)
178       return;
179  
180     while((blres = sroi_getnext_build(fname)) > 0)
181     {
182         input_model(fname);
183         if ((model = get_model()) == NULL)
184             return;
185         
186         m = model->m->el.double_v;
187         p = model->profile->el.float_v;
188         pg = model->g_profile->el.double_v;
189 
190         if (i == 0)
191         {
192             length = 2 * model->vsize - 1;
193             plength = 2 * model->pcavsize - 1;
194             start = model->vsize - model->pcavsize;
195             NX = model->c;
196             NXY = 2 * NX;
197             ng = NX * length;
198             M_mean = vector_alloc(NXY, double_v);
199             M_cov = matrix_alloc(NXY, NXY, matrix_full, double_v);
200             M_Evec = matrix_alloc(NXY, NXY, matrix_full, double_v);
201             M_Eval = vector_alloc(NXY, double_v);
202             M_xy = vector_alloc(NXY, double_v);
203             mm = (double *) M_mean->data;
204             mc = M_cov->el.double_v;
205             mxy = (double *) M_xy->data;
206 
207             G_mean = vector_alloc(ng, double_v);
208             G_cov = matrix_alloc(ng, ng, matrix_full, double_v);
209             G_Evec = matrix_alloc(ng, ng, matrix_full, double_v);
210             G_Eval = vector_alloc(ng, double_v);
211             R_Evec = matrix_alloc(ng, NX, matrix_full, double_v);
212             R_Eval = vector_alloc(NX, double_v);
213             gm = (double *) G_mean->data;
214             gc = G_cov->el.double_v;
215 
216             P_mean = (Vector **) pvector_alloc(0, model->c);
217             P_cov = (Matrix **) pvector_alloc(0, model->c);
218             pm = (double **) pvector_alloc(0, model->c);
219             pc = (double ***) pvector_alloc(0, model->c);
220             for (r = 0; r < model->c; r++)
221             {
222                 P_mean[r] = vector_alloc(length, double_v);
223                 pm[r] = (double *) P_mean[r]->data;
224                 P_cov[r] = matrix_alloc(plength, plength,
225                                         matrix_full, double_v);
226                 pc[r] = P_cov[r]->el.double_v;
227             }
228             P_Evec = matrix_alloc(plength, plength,
229                                   matrix_full, double_v);
230             P_Eval = vector_alloc(plength, double_v);
231         }
232         
233         for (r = 0; r < NX; r++)
234         {
235             mxy[r] = m[3][r];   
236             mxy[NX + r] = m[2][r];
237         }
238         for (r = 0; r < NXY; r++)
239             mm[r] += mxy[r];
240         for (r = 0; r < NXY; r++)
241         {
242             mc[r][r] += mxy[r] * mxy[r];
243             for (c = r + 1; c < NXY; c++)
244                 mc[r][c] += mxy[r] * mxy[c];
245         }
246 
247         if (gradient_on)
248         {
249             for (r = 0; r < model->c; r++)
250             {
251                 pg[r][0] = (double) (p[r][1] - p[r][0]);
252                 pg[r][1] = (double) (p[r][2] - p[r][0]);
253                 for (c = 2; c < length - 2; c++)
254                     pg[r][c] = (double) (p[r][c + 2] - p[r][c - 2]);
255                 pg[r][length - 2] = (double) (p[r][length - 1] - p[r][length - 3]);
256                 pg[r][length - 1] = (double) (p[r][length - 1] - p[r][length - 2]);
257             }
258         } else
259         {
260             for (r = 0; r < model->c; r++)
261                 for (c = 0; c < length; c++)
262                     pg[r][c] = (double) p[r][c];
263         }
264         
265         if (normalise_on)
266         {
267             for (j = 0; j < model->c; j++)
268             {
269                 for (r = 0, norm = 0; r < length; r++)
270                     norm += fabs(pg[j][r]);
271                 norm /= (double) length;
272                 for (r = 0; r < length; r++)
273                     pg[j][r] /= norm;
274             }
275         }
276         
277         if (global_on)
278         {
279             for (r = 0, j = 0; r < model->c; r++)
280                 for (c = 0; c < length; c++)
281                     gm[j++] += pg[r][c];
282             for (r = 0; r < ng; r++)
283             {
284                 gc[r][r] += gm[r] * gm[r];
285                 for (c = r + 1; c < ng; c++)
286                     gc[r][c] += gm[r] * gm[c];
287             }
288         } else
289         {
290             for (j = 0; j < model->c; j++)
291                 for (r = 0; r < length; r++)
292                     pm[j][r] += pg[j][r];
293             for (j = 0; j < model->c; j++)
294             {
295                 for (r = 0; r < plength; r++)
296                 {
297                     pc[j][r][r] += pg[j][start + r] * pg[j][start + r];
298                     for (c = r + 1; c < plength; c++)
299                         pc[j][r][c] += pg[j][start + r] * pg[j][start + c];
300                 }
301             }
302         }
303         tx += model->tx;
304         ty += model->ty;
305         s += model->s;
306         theta += model->theta;
307         printf("%d %s\n", i++, fname);
308     }
309 
310     sroi_close_buildlist();
311    
312     if (blres != 0)
313     {
314       error("pca_model: buildlist file not finished - aborting", warning);
315       return;
316     }
317   
318     for (r = 0; r < NXY; r++)
319         mm[r] /= n;
320     for (r = 0; r < NXY; r++)
321     {
322         for (c = r; c < NXY; c++)
323             mc[r][c] = mc[r][c] / n - mm[r] * mm[c];
324         for (c = r + 1; c < NXY; c++)
325             mc[c][r] = mc[r][c];
326     }
327     for (r = 0; r < NX; r++)
328     {
329         m[3][r] = mm[r];
330         m[2][r] = mm[NX + r];
331     }
332     
333     if (global_on)
334     {
335         for (r = 0; r < ng; r++)
336             gm[r] /= n;
337         for (r = 0; r < ng; r++)
338         {
339             for (c = r; c < ng; c++)
340                 gc[r][c] = gc[r][c] / n - gm[r] * gm[c];
341             for (c = r + 1; c < ng; c++)
342                 gc[c][r] = gc[r][c];
343         }
344         for (r = 0, j = 0; r < NX; r++)
345         {
346             for (c = 0; c < length; c++)
347                 p[r][c] = (int) gm[j++];
348         }
349     } else
350     {
351         for (j = 0; j < model->c; j++)
352             for (r = 0; r < length; r++)
353                 pm[j][r] /= n;
354         for (j = 0; j < model->c; j++)
355         {
356             for (r = 0; r < plength; r++)
357             {
358                 for (c = r; c < plength; c++)
359                     pc[j][r][c] = pc[j][r][c] / n - pm[j][start + r] * pm[j][start + c];
360                 for (c = r + 1; c < plength; c++)
361                     pc[j][c][r] = pc[j][r][c];
362             }
363         }
364         for (r = 0; r < model->c; r++)
365         {
366             for (c = 0; c < length; c++)
367                 p[r][c] = (float) pm[r][c];
368         }
369     }
370 
371     model->tx = tx / n;
372     model->ty = ty / n;
373     model->s = s / n;
374     model->theta = theta / n;
375 
376     output_model(meanfn);
377 
378     if ((f = fopen(pcafn, "wt")) == NULL)
379     {
380         string_append(temps, "can't open file", pcafn, 0);
381         error(temps, non_fatal);
382         return;
383     }
384     matrix_eigen_sym(M_cov, M_Eval, M_Evec);
385 
386     matrix_swap_all_cols(M_Evec);
387     vector_swap(M_Eval);
388     dvector_sqrt(M_Eval);
389     vector_format(M_Eval);
390     dmatrix_print(f, M_Evec);
391     dvector_printf(f, M_Eval);
392     if (global_on)
393     {
394         matrix_eigen_sym(G_cov, G_Eval, G_Evec);
395         matrix_swap_all_cols(G_Evec);
396         vector_swap(G_Eval);
397         dvector_sqrt(G_Eval);
398         vector_format(G_Eval);
399         for (c = 0; c < model->c; c++)
400         {
401             for (r = 0; r < ng; r++)
402                 R_Evec->el.double_v[r][c] = G_Evec->el.double_v[r][c];
403             ((double *) (R_Eval->data))[c] = ((double *) (G_Eval->data))[c];
404         }
405         dmatrix_print(f, R_Evec);
406         dvector_printf(f, R_Eval);
407     } else
408     {
409         for (j = 0; j < model->c; j++)
410         {
411             matrix_eigen_sym(P_cov[j], P_Eval, P_Evec);
412             matrix_swap_all_cols(P_Evec);
413             vector_swap(P_Eval);
414             dvector_sqrt(P_Eval);
415             dmatrix_print(f, P_Evec);
416             dvector_printf(f, P_Eval);
417         }
418     }
419     vector_free(M_mean);
420     matrix_free(M_cov);
421     vector_free(M_Eval);
422     matrix_free(M_Evec);
423     vector_free(G_mean);
424     matrix_free(G_cov);
425     vector_free(G_Eval);
426     matrix_free(G_Evec);
427     vector_free(R_Eval);
428     matrix_free(R_Evec);
429     vector_free(P_Eval);
430     matrix_free(P_Evec);
431     for (i = 0; i < model->c; i++)
432     {
433         vector_free(P_mean[i]);
434         matrix_free(P_cov[i]);
435     }
436     fclose(f);
437 }
438 

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