1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgEM_mix.c,v $
37 * Date : $Date: 2005/07/03 23:46:54 $
38 * Version : $Revision: 1.10 $
39 * CVS Id : $Id: imgEM_mix.c,v 1.10 2005/07/03 23:46:54 paul Exp $
40 */
41
42 /**
43 * @file
44 * @brief EM Algorithm: Gaussian distribution model with linear partial volume terms.
45 *
46 */
47
48 #include "imgEM_mix.h"
49
50 #if HAVE_CONFIG_H
51 #include <config.h>
52 #endif
53
54
55 #include <stdio.h>
56 #include <math.h>
57 #include <tina/sys/sysDef.h>
58 #include <tina/sys/sysPro.h>
59 #include <tina/math/mathDef.h>
60 #include <tina/math/mathPro.h>
61 #include <tina/file/filePro.h>
62 #include <tina/image/imgDef.h>
63 #include <tina/image/img_EMDef.h>
64
65 #ifdef _PCC
66 #ifndef MINGW
67 #include <tina/sys/syswin_support.h> // for erf()
68 #endif
69 #endif
70
71 #define GMAX2 49.0
72
73
74 static double unit_gaussian2(double z)
75 {
76 double x;
77 if (z > GMAX2)
78 return (0.0);
79 x = exp(-z / 2.0);
80 return (x);
81 }
82
83
84 Mixmodel *mixmodel_alloc(int ndim, int nmix, double background)
85 {
86 Mixmodel *model = NULL;
87 int i, k;
88
89 model = (Mixmodel *) ralloc(sizeof(Mixmodel));
90 model->nmix = nmix;
91 model->ndim = ndim;
92 model->background = background;
93 model->vectors = (double **) pvector_alloc(0, model->nmix);
94 model->normal = (double *) nvector_alloc(0, model->nmix, sizeof(double));
95 model->volume = (double *) nvector_alloc(0, model->nmix, sizeof(double));
96 model->gradscale = NULL;
97 for (i = 0; i < model->nmix; i++)
98 {
99 model->vectors[i] = (double *) nvector_alloc(0, model->ndim, sizeof(double));
100 }
101 model->cov = (double ***) pvector_alloc(0, model->nmix);
102 for (k = 0; k < model->nmix; k++)
103 {
104 model->cov[k] = darray_alloc(0, 0, model->ndim, model->ndim);
105 }
106 model->par_dens = darray_alloc(0, 0, model->nmix, model->nmix);
107
108 model->tissue_type = (char **) pvector_alloc(0, model->nmix);
109 for (i = 0; i < model->nmix; i++)
110 {
111 model->tissue_type[i] = (char *) nvector_alloc(0, 24, sizeof(char));
112 }
113 return (model);
114 }
115
116 void mixmodel_free(Mixmodel * model)
117 {
118 int i, k;
119 for (i = 0; i < model->nmix; i++)
120 {
121 nvector_free(model->tissue_type[i], 0, sizeof(char));
122 }
123 pvector_free(model->tissue_type, 0);
124 darray_free(model->par_dens, 0, 0, model->nmix, model->nmix);
125 for (k = 0; k < model->nmix; k++)
126 {
127 darray_free(model->cov[k], 0, 0, model->ndim, model->ndim);
128 }
129 pvector_free(model->cov, 0);
130 for (i = 0; i < model->nmix; i++)
131 {
132 nvector_free(model->vectors[i], 0, sizeof(double));
133 }
134 pvector_free(model->vectors, 0);
135 nvector_free(model->normal, 0, sizeof(double));
136 nvector_free(model->volume, 0, sizeof(double));
137 if (model->gradscale != NULL)
138 model->gradscale = darray_free(model->gradscale,0,0,model->nmix,model->nmix);
139 rfree(model);
140 }
141
142 Mixmodel *input_mixmodel(char *file_name)
143 {
144 Mixmodel *model = NULL;
145 FILE *stream = NULL;
146 int i, j, k, nmix, ndim;
147 double c, grad;
148 double background;
149
150
151 stream = (FILE *) fopen_2(file_name, "r");
152 if (stream == NULL)
153 {
154 error("Non-existant parameter file", warning);
155 return (NULL);
156 }
157 /*** nmix - number of mixtures (classes) ***/
158 /*** ndim - number of dimensions (images) ***/
159 /*** background - term for background (noise) ***/
160
161 fscanf(stream, "%d %d %lf", &nmix,
162 &ndim,
163 &background);
164
165 if (model != NULL)
166 mixmodel_free(model);
167 model = mixmodel_alloc(ndim, nmix, background);
168
169 for (i = 0; i < model->nmix; i++)
170 {
171 for (j = 0; j < model->ndim; j++)
172 {
173 fscanf(stream, " %lf", &(model->vectors[i][j]));
174 }
175 fscanf(stream, " %lf ", &(model->volume[i]));
176 }
177
178 for (k = 0; k < model->nmix; k++)
179 {
180 model->volume[k] = 1.0;
181 for (i = 0; i < model->ndim; i++)
182 {
183 for (j = 0; j < model->ndim; j++)
184 {
185 fscanf(stream, " %lf", &(model->cov[k][i][j]));
186 model->cov[k][i][j] *= fabs(model->cov[k][i][j]);
187 }
188 model->volume[k] *= sqrt(TWOPI);
189 }
190 c = darray_det(model->cov[k], model->ndim);
191 model->volume[k] /= sqrt(c);
192 }
193
194 for (i = 0; i < model->nmix; i++)
195 {
196 for (j = 0; j < model->nmix; j++)
197 {
198 fscanf(stream, " %lf", &(model->par_dens[i][j]));
199 if (i == j)
200 model->normal[i] = model->par_dens[i][j];
201 }
202 }
203 format("\ntissue types = ");
204 for (i = 0; i < model->nmix; i++)
205 {
206 fscanf(stream, " %s", model->tissue_type[i]);
207 format("%s ", model->tissue_type[i]);
208 }
209 format("\n",0);
210 if (fscanf(stream, " %lf", &grad) != EOF)
211 {
212 model->gradscale = darray_alloc(0,0,model->nmix,model->nmix);
213 for (i = 0; i < model->nmix; i++)
214 {
215 for (j = 0; j < model->nmix; j++)
216 {
217 model->gradscale[i][j]= grad;
218 fscanf(stream, " %lf", &grad);
219 }
220 }
221 }
222 fclose(stream);
223 return (model);
224 }
225
226 void output_mixmodel(Mixmodel * model, char *file_name)
227 {
228 FILE *stream = NULL;
229 double tmp;
230 int i, j, k;
231
232 stream = (FILE *) fopen_2(file_name, "w");
233 if (stream == NULL)
234 return;
235
236 fprintf(stream, "%d %d %18.14f \n", model->nmix, model->ndim, model->background);
237
238 for (i = 0; i < model->nmix; i++)
239 {
240 for (j = 0; j < model->ndim; j++)
241 {
242 fprintf(stream, " %f", model->vectors[i][j]);
243 }
244 fprintf(stream, " %f \n", model->volume[i]);
245 }
246
247 for (k = 0; k < model->nmix; k++)
248 {
249 for (i = 0; i < model->ndim; i++)
250 {
251 for (j = 0; j < model->ndim; j++)
252 {
253 tmp = model->cov[k][i][j];
254 if (tmp != 0.0)
255 tmp /= sqrt(fabs(tmp));
256 fprintf(stream, " %f", tmp);
257 }
258 fprintf(stream, " \n");
259 }
260 }
261
262 for (i = 0; i < model->nmix; i++)
263 {
264 for (j = 0; j < model->nmix; j++)
265 {
266 fprintf(stream, " %f", model->par_dens[i][j]);
267 }
268 fprintf(stream, " \n");
269 }
270
271 for (i = 0; i < model->nmix; i++)
272 fprintf(stream, " %s", model->tissue_type[i]);
273 fprintf(stream, " \n");
274
275 if (model->gradscale != NULL)
276 {
277 for (i = 0; i < model->nmix; i++)
278 {
279 for (j = 0; j < model->nmix; j++)
280 {
281 fprintf(stream, " %f", model->gradscale[i][j]);
282 }
283 fprintf(stream, " \n");
284 }
285 }
286
287
288 fclose(stream);
289 }
290
291
292 Mixmodel *get_submodel(Mixmodel * model, int n1, int n2)
293 {
294 int i, j, k, p, r, n;
295 Mixmodel *model_2d = NULL;
296
297 if (model == NULL)
298 {
299 error("Input file is missing!", warning);
300 return (NULL);
301 }
302 model_2d = mixmodel_alloc(n2 - n1 + 1, model->nmix, model->background);
303
304 for (i = 0; i < model->nmix; i++)
305 {
306 n = n1;
307 for (j = 0; j < model_2d->ndim; j++)
308 {
309 model_2d->vectors[i][j] = model->vectors[i][n];
310 n++;
311 }
312 model_2d->volume[i] = model->volume[i];
313 }
314 for (k = 0; k < model->nmix; k++)
315 {
316 model_2d->volume[k] = 1.0;
317 p = n1;
318 for (i = 0; i < model_2d->ndim; i++)
319 {
320 r = n1;
321 for (j = 0; j < model_2d->ndim; j++)
322 {
323 model_2d->cov[k][i][j] = model->cov[k][p][r];
324 r++;
325 }
326 p++;
327 model_2d->volume[k] *= sqrt(TWOPI);
328 }
329 model_2d->volume[k] /= sqrt(darray_det(model_2d->cov[k], model_2d->ndim));
330 }
331
332 for (i = 0; i < model->nmix; i++)
333 {
334 for (j = 0; j < model->nmix; j++)
335 {
336 model_2d->par_dens[i][j] = model->par_dens[i][j];
337 if (i == j)
338 model_2d->normal[i] = model->par_dens[i][j];
339 }
340 }
341 for (i = 0; i < model_2d->nmix; i++)
342 {
343 strcpy(model_2d->tissue_type[i], model->tissue_type[i]);
344 }
345 return (model_2d);
346 }
347
348 double pure_density(Mixmodel * m, float *x, int blob)
349 {
350
351 int i, j;
352 double z, tmp, density = 1.0;
353 double **covar, *vblob;
354
355 z = 0.0;
356 covar = m->cov[blob];
357 vblob = m->vectors[blob];
358 for (j = 0; j < m->ndim; j++)
359 {
360 tmp = 0.0;
361 for (i = 0; i < m->ndim; i++)
362 {
363 tmp += (x[i] - vblob[i]) * covar[i][j];
364 }
365 z += tmp * (x[j] - vblob[j]);
366 }
367
368 density = unit_gaussian2(z) / m->volume[blob];
369 return (density);
370 }
371
372 void get_covar(Mixmodel * m, int blob1, int blob2, double a, double **newcov)
373 {
374 int i, j;
375 double **cov1, **cov2;
376 cov1 = m->cov[blob1];
377 cov2 = m->cov[blob2];
378 for (j = 0; j < m->ndim; j++)
379 {
380 for (i = 0; i < m->ndim; i++)
381 {
382 newcov[j][i] = a * cov1[j][i]
383 + (1.0 - a) * cov2[j][i];
384 }
385 }
386 }
387
388 double linear_fraction(Mixmodel * m, float *v, int blob1, int blob2, float a)
389 {
390 int i, j;
391 double xs = 0.0, x = 0.0;
392 static double **newcov = NULL;
393 static int dim = -1;
394 double vec1, vec2, *vblob1, *vblob2;
395 double result;
396
397 if (dim != m->ndim && newcov)
398 darray_free(newcov, 0, 0, dim, dim);
399 dim = m->ndim;
400 if (!newcov)
401 newcov = darray_alloc(0, 0, dim, dim);
402
403 get_covar(m, blob1, blob2, a, newcov);
404
405 vblob1 = m->vectors[blob1];
406 vblob2 = m->vectors[blob2];
407 for (j = 0; j < dim; j++)
408 {
409 vec1 = 0.0;
410 vec2 = 0.0;
411 for (i = 0; i < dim; i++)
412 {
413 vec1 += (vblob1[i] - vblob2[i]) *
414 newcov[i][j];
415 vec2 += (v[i] - vblob2[i]) * newcov[i][j];
416 }
417 xs += vec1 * (vblob1[j] - vblob2[j]);
418 x += vec2 * (vblob1[j] - vblob2[j]);
419 }
420 result = x / xs; /* normal distance from the mean of blob2 */
421
422 if (result < 0.0)
423 {
424 result = 0.0;
425 }
426 if (result > 1.0)
427 {
428 result = 1.0;
429 }
430 if (result > 0.0 || result <= 0.0); /* NaN trap */
431 else
432 result = 0.5;
433 return (result);
434 }
435
436 double conv_gauss_tri(double sig, double xs, double xl, double x)
437 {
438 double xdiff, a, b, c, m, lin;
439 double xa, xb, ans1, ans2, ans;
440
441 a = xs;
442 b = xl;
443 if (xs > xl)
444 {
445 b = xs;
446 a = xl;
447 }
448 /*
449 h = 1.0/(b-a);
450 */
451 xdiff = xl - xs;
452 m = 1.0 / xdiff;
453 c = -xs * m;
454
455 xa = (x - a) / (sqrt(2.0) * sig);
456 xb = (x - b) / (sqrt(2.0) * sig);
457 lin = m * x + c;
458
459 ans1 = -0.5 * lin * (erf(xb) - erf(xa));
460 ans2 = exp(-xa * xa) - exp(-xb * xb);
461 ans2 *= m * sig / sqrt(TWOPI);
462 ans = ans1 + ans2;
463 return ans;
464 }
465
466 double part_fraction(Mixmodel * m, float *v, int blob1, int blob2)
467 {
468 double a;
469
470 a = linear_fraction(m, v, blob1, blob2, 0.5);
471 a = linear_fraction(m, v, blob1, blob2, (float) a);
472 a = linear_fraction(m, v, blob1, blob2, (float) a);
473 a = linear_fraction(m, v, blob1, blob2, (float) a);
474 a = linear_fraction(m, v, blob1, blob2, (float) a);
475 a = linear_fraction(m, v, blob1, blob2, (float) a);
476
477 return(a);
478 }
479
480 double part_density(Mixmodel * m, float *v, int blob1, int blob2, double a)
481 {
482 int i, j;
483 double density;
484 static double **newcov = NULL;
485 static int dim = -1;
486 double newnorm;
487 double y1, x1 = 0.0, n1 = 0.0, r1 = 0.0;
488 double vec1, vec2, *vblob1, *vblob2;
489
490 if (blob1 == blob2 || m == NULL)
491 return (0.0);
492 vblob1 = m->vectors[blob1];
493 vblob2 = m->vectors[blob2];
494
495 if (dim != m->ndim && newcov)
496 darray_free(newcov, 0, 0, dim, dim);
497 dim = m->ndim;
498 if (!newcov)
499 newcov = darray_alloc(0, 0, dim, dim);
500
501 get_covar(m, blob1, blob2, a, newcov);
502 newnorm = exp(dim * log(TWOPI) / 2.0) / sqrt(darray_det(newcov, dim));
503
504 for (j = 0; j < dim; j++)
505 {
506 vec1 = 0.0;
507 vec2 = 0.0;
508 for (i = 0; i < dim; i++)
509 {
510 vec1 += (vblob1[i] - vblob2[i]) * newcov[i][j];
511 vec2 += (v[i] - vblob1[i]) * newcov[i][j];
512 }
513 r1 += vec2 * (v[j] - vblob1[j]);
514 n1 += vec1 * (vblob1[j] - vblob2[j]);
515 x1 += vec2 * (vblob2[j] - vblob1[j]);
516 }
517 n1 = sqrt(n1);
518 x1 = x1 / n1;
519 if (x1 > 0.0 || x1 <= 0.0); /* NaN trap */
520 else
521 {
522 x1 = 0.0;
523 n1 = 0.0;
524 error("Class overlap", warning);
525 }
526 y1 = (r1 - x1 * x1);
527 if (y1 < 0.0)
528 y1 = 0.0;
529
530 density = unit_gaussian2(y1);
531 /*
532 density *= 2.0*nmr_gaus_tri(1.0,n1,0.0,x1)*sqrt(TWOPI)/(n1*newnorm);
533 */
534 density *= 2.0 * conv_gauss_tri(1.0, n1, 0.0, x1) * sqrt(TWOPI) / (n1 * newnorm);
535 darray_free(newcov, 0, 0, dim, dim);
536 if(density <= 0.0 ) density =0.0;
537 return (density);
538 }
539
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.