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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.