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

Linux Cross Reference
Tina5/tina-libs/tina/image/imgEM_mix.c

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

  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 

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