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

Linux Cross Reference
Tina5/tina-libs/tina/image/imgEM_max.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_max.c,v $
 37  * Date    :  $Date: 2003/09/22 16:09:02 $
 38  * Version :  $Revision: 1.7 $
 39  * CVS Id  :  $Id: imgEM_max.c,v 1.7 2003/09/22 16:09:02 tony Exp $
 40  */
 41 
 42 /** 
 43  *  @file
 44  *  @brief EM Algorithm: Maximisation (M) step
 45  *
 46  * Maximisation (M) step of EM algorithm     
 47  * Calculation of model parameters ; mean values covariances, pure and
 48  * partial volume proportions.
 49  *
 50  */
 51 
 52 #include "imgEM_max.h"
 53 
 54 #if HAVE_CONFIG_H
 55 #include <config.h>
 56 #endif
 57 
 58 
 59 #include <stdio.h>
 60 #include <math.h>
 61 #include <tina/sys/sysDef.h>
 62 #include <tina/sys/sysPro.h>
 63 #include <tina/math/mathDef.h>
 64 #include <tina/math/mathPro.h>
 65 #include <tina/image/imgDef.h>
 66 #include <tina/image/imgPro.h>
 67 #include <tina/image/img_EMDef.h>
 68 #include <tina/image/imgEM_probs.h>
 69 #include <tina/image/imgEM_estep.h>
 70 
 71 
 72 void         ***seq_limits(int *lz, int *uz, int *ly, int *uy, int *lx, int *ux);
 73 
 74 /*** Calculation of inverse of covariance matrix ***/
 75 
 76 void            cov_update(void ***pim_array, Imregion * roi, Mixmodel * m)
 77 {
 78         int             i, j, k, loop1, loop2, loop3;
 79         double          tot_pix, this_pix, ***C, ***Ctot;
 80         float          *data, *row;
 81         void         ***imptrs = NULL;
 82         int             imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
 83         int             lx, ux, ly, uy;
 84         Imrect         *im_tot = em_get_im_tot();
 85 
 86         if (pim_array == NULL)
 87         {
 88                 error("Recalculate densities - <E-step>", warning);
 89                 return;
 90         }
 91         imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
 92                             &imptrux);
 93         data = (float *) fvector_alloc(imptrlz, imptruz);
 94 
 95         C = (double ***) pvector_alloc(0, m->nmix);
 96         for (k = 0; k < m->nmix; k++)
 97                 C[k] = darray_alloc(0, 0, m->ndim, m->ndim);
 98 
 99         Ctot = (double ***) pvector_alloc(0, m->nmix);
100         for (k = 0; k < m->nmix; k++)
101                 Ctot[k] = darray_alloc(0, 0, m->ndim, m->ndim);
102 
103         for (i = 0; i < m->nmix; i++)
104                 for (j = 0; j < m->ndim; j++)
105                         for (k = 0; k < m->ndim; k++)
106                         {
107                                 C[i][j][k] = 0.0;
108                                 Ctot[i][j][k] = 0.0;
109                         }
110         lx = roi->lx;
111         ux = roi->ux;
112         ly = roi->ly;
113         uy = roi->uy;
114 
115         for (i = ly; i < uy; ++i)
116         {
117                 for (j = lx; j < ux; ++j)
118                 {
119 
120                         tot_pix = im_get_pixf(im_tot, i, j);
121                         for (k = imptrlz; k < imptruz; k++)
122                         {
123                                 row = imptrs[k][i];
124                                 data[k] = row[j];
125                         }
126                         /* Covariance calculation */
127                         for (loop1 = 0; loop1 < m->nmix; loop1++)
128                         {
129                                 this_pix = im_get_pixf((Imrect *) pim_array[loop1][loop1], i, j);
130 
131                                 for (loop2 = 0; loop2 < m->ndim; loop2++)
132                                 {
133                                         for (loop3 = 0; loop3 < m->ndim; loop3++)
134                                         {
135                                                 if (im_get_pixf(im_tot, i, j) > 0.0)
136                                                 {
137                                                         C[loop1][loop2][loop3] += (data[loop2 + imptrlz] - m->vectors[loop1][loop2])
138                                                                 * (data[loop3 + imptrlz] - m->vectors[loop1][loop3])
139                                                                 * this_pix / tot_pix;
140 
141                                                         Ctot[loop1][loop2][loop3] += this_pix / tot_pix;
142                                                 }
143                                         }
144                                 }
145                         }
146 
147                 }               /* end of j loop */
148         }                       /* end of i loop */
149 
150         for (k = 0; k < m->nmix; k++)
151         {
152                 for (i = 0; i < m->ndim; i++)
153                 {
154                         for (j = 0; j < m->ndim; j++)
155                                 if (Ctot[k][i][j] > 0.0)
156                                         C[k][i][j] = C[k][i][j] / Ctot[k][i][j];
157                 }
158 
159                 if (darray_invert(C[k], m->ndim) == NULL)
160                 {
161                         error("Not enough data for number of classes specified!", warning);
162                         return;
163                 } else
164                         m->cov[k] = darray_invert(C[k], m->ndim);
165 
166                 m->volume[k] = exp(m->ndim * log(TWOPI) / 2.0) / sqrt(darray_det(m->cov[k], m->ndim));
167         }
168 
169         /*
170            my_log_proc();
171         */
172         for (k = 0; k < m->nmix; k++)
173         {
174                 darray_free(C[k], 0, 0, m->ndim, m->ndim);
175                 darray_free(Ctot[k], 0, 0, m->ndim, m->ndim);
176         }
177         pvector_free(C, 0);
178         pvector_free(Ctot, 0);
179 }
180 
181 /*** Caclulation of priors ***/
182 void            priors_update(void ***pim_array, Imregion * roi, Mixmodel * m)
183 {
184         int             i, j, loop, loop1, t = 0;
185         double        **f;
186         int             lx, ux, ly, uy;
187         Imrect         *im_tot = em_get_im_tot();
188 
189         if (pim_array == NULL)
190         {
191                 error("Recalculate densities - <E-step>", warning);
192                 return;
193         }
194         f = darray_alloc(0, 0, m->nmix, m->nmix);
195         for (i = 0; i < m->nmix; i++)
196                 for (j = 0; j < m->nmix; j++)
197                         f[i][j] = 0.0;
198         lx = roi->lx;
199         ux = roi->ux;
200         ly = roi->ly;
201         uy = roi->uy;
202         for (loop = 0; loop < m->nmix; loop++)
203         {
204                 for (loop1 = 0; loop1 < m->nmix; loop1++)
205                 {
206                         if (pim_array[loop][loop1] != NULL)
207                         {
208                                 for (i = ly; i < uy; ++i)
209                                 {
210                                         for (j = lx; j < ux; ++j)
211                                         {
212                                                 /*
213                                                  * Calculate volume
214                                                  * normalization factors
215                                                  * (priors)
216                                                  */
217 
218                                                 if (im_get_pixf(im_tot, i, j) > 0.0)
219                                                         f[loop][loop1]
220                                                                 += im_get_pixf((Imrect *) pim_array[loop][loop1], i, j)
221                                                                 / im_get_pixf(im_tot, i, j);
222                                                 else
223                                                         t += 1;
224                                         }       /* end of j loop */
225                                 }       /* end of i loop */
226                         }
227                 }
228         }
229         if (t > 0)
230                 format("\n no of pixels not labelled = %i", t);
231 
232         for (i = 0; i < m->nmix; ++i)
233         {
234                 for (j = 0; j < m->nmix; ++j)
235                 {
236                         if (i == j)
237                         {
238                                 m->par_dens[i][j] = f[i][j];
239                                 m->normal[i] = f[i][i];
240                         } else
241                                 m->par_dens[i][j] = (f[i][j] + f[j][i]) * 0.5;
242                 }
243         }
244         /*
245            my_log_proc();
246         */
247         darray_free(f, 0, 0, m->nmix, m->nmix);
248 }
249 
250 /*** Calculation of means ***/
251 
252 void            mean_update(void ***pim_array, Imregion * roi, Mixmodel * m)
253 {
254         Imrect         *im = NULL;
255         int             i, j, k, loop, loop1;
256         double          tot, this_pix, **M, **Mtot;
257         float          *data, *row;
258         void         ***imptrs = NULL;
259         int             imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
260         int             lx, ux, ly, uy;
261         Imrect         *im_tot = em_get_im_tot();
262 
263         if (im_tot == NULL)
264                 im_tot = prob_im_tot(m, pim_array, roi);
265 
266         if (pim_array == NULL)
267         {
268                 error("Recalculate densities - <E-step>", warning);
269                 return;
270         }
271         imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
272                             &imptrux);
273         data = (float *) fvector_alloc(imptrlz, imptruz);
274 
275         lx = roi->lx;
276         ux = roi->ux;
277         ly = roi->ly;
278         uy = roi->uy;
279 
280         M = darray_alloc(0, 0, m->nmix, m->ndim);
281         Mtot = darray_alloc(0, 0, m->nmix, m->ndim);
282         for (i = 0; i < m->nmix; i++)
283                 for (j = 0; j < m->ndim; j++)
284                 {
285                         M[i][j] = 0.0;
286                         Mtot[i][j] = 0.0;
287                 }
288         for (i = ly; i < uy; ++i)
289         {
290                 for (j = lx; j < ux; ++j)
291                 {
292                         for (k = imptrlz; k < imptruz; k++)
293                         {
294                                 row = imptrs[k][i];
295                                 data[k] = row[j];
296                         }
297 
298                         /*
299                          * Calculate mean value for each pixel in each
300                          * density image
301                          */
302 
303                         tot = im_get_pixf(im_tot, i, j);
304                         for (loop = 0; loop < m->nmix; loop++)
305                                 for (loop1 = 0; loop1 < m->ndim; loop1++)
306                                 {
307                                         im = (Imrect *) pim_array[loop][loop];
308                                         this_pix = im_get_pixf(im, i, j);
309                                         if (tot > 0.0)
310                                         {
311                                                 M[loop][loop1] += this_pix * data[loop1 + imptrlz] / tot;
312                                                 Mtot[loop][loop1] += this_pix / tot;
313                                         }
314                                 }
315                 }               /* end of j loop */
316         }                       /* end of i loop */
317         for (i = 0; i < m->nmix; i++)
318                 for (j = 0; j < m->ndim; j++)
319                         if (Mtot[i][j] > 0.0)
320                                 m->vectors[i][j] = M[i][j] / Mtot[i][j];
321         /*
322            my_log_proc();
323         */
324         darray_free(M, 0, 0, m->nmix, m->ndim);
325         darray_free(Mtot, 0, 0, m->nmix, m->ndim);
326 }
327 

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