1 /**********
2 *
3 * This file is part of the TINA Open Source Image Analysis Environment
4 * henceforth known as TINA
5 *
6 * TINA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as
8 * published by the Free Software Foundation.
9 *
10 * TINA is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with TINA; if not, write to the Free Software Foundation, Inc.,
17 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 **********
20 *
21 * Program : TINA
22 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgEM_probs.c,v $
23 * Date : $Date: 2005/01/23 19:10:21 $
24 * Version : $Revision: 1.10 $
25 * CVS Id : $Id: imgEM_probs.c,v 1.10 2005/01/23 19:10:21 paul Exp $
26 *
27 * Notes :
28 *
29 * Gio Buonaccorsi commented out prob_all as it is required only for testing.
30 * 19 Feb 2003.
31 *
32 **********************************************
33 * em_params_max.c *
34 **********************************************
35 * Maximisation (M) step of EM algorithm *
36 * Calculation of model parameters *
37 **********************************************
38 *
39 *********
40 */
41
42 #include "imgEM_probs.h"
43 #include "imgEM_grad.h"
44
45 #if HAVE_CONFIG_H
46 #include <config.h>
47 #endif
48
49 #include <stdio.h>
50 #include <math.h>
51 #include <string.h>
52 #include <tina/sys/sysDef.h>
53 #include <tina/sys/sysPro.h>
54 #include <tina/math/mathDef.h>
55 #include <tina/math/mathPro.h>
56 #include <tina/image/imgDef.h>
57 #include <tina/image/imgPro.h>
58 //#include "imgEM_mix.h"
59 //#include "img_EMDef.h"
60
61
62 extern List *get_current_seq_current_el(void);
63 extern void ***seq_limits(int *lz, int *uz, int *ly, int *uy, int *lx, int *ux);
64 extern double **get_kgrad_arr(void);
65
66 Imrect *prob_im_grad( Mixmodel *m, Imrect *im_prob, Imregion *roi, int select1, int select2, Imrect *im_slope)
67 {
68 int lx, ux, ly, uy;
69 int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
70 double tissue;
71 float *row, *row1, *data;
72 int i,j,k,loop,loop1;
73 void ***imptrs = NULL;
74 double **Sk = NULL, y, a, d;
75 double s0=1.0,Pgrad=0.0,sf=0.0;
76 double fract1, fract2;
77
78 if (im_prob == NULL)
79 {
80 error("Density image memory not allocated!", warning);
81 return(NULL);
82 }
83 if (im_slope == NULL)
84 {
85 error("Gradient image not calculated!", warning);
86 return(NULL);
87 }
88
89 if (m == NULL)
90 {
91 error("Input file is missing!", warning);
92 return(NULL);
93 }
94 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly,
95 &imptruy, &imptrlx, &imptrux);
96 data = (float *)fvector_alloc(imptrlz, imptruz);
97
98
99 if (imptrlx > (lx = roi->lx)) lx = imptrlx;
100 if (imptrux < (ux = roi->ux)) ux = imptrux;
101 if (imptrly > (ly = roi->ly)) ly = imptrly;
102 if (imptruy < (uy = roi->uy)) uy = imptruy;
103
104
105
106 s0 = (double)sqrt(m->ndim);
107 /* Get the values of scaling factors Sk */
108 Sk = m->gradscale;
109 if(Sk == NULL)
110 {
111 Sk = darray_alloc(0,0,m->nmix,m->nmix);
112 m->gradscale = Sk;
113 for(loop = 0; loop < m->nmix; loop++)
114 for(loop1 = 0; loop1 < m->nmix; loop1++)
115 Sk[loop][loop1] = 1.0;
116 }
117 /* Calculate normalisation factors for gradient distribution */
118
119 row1 = fvector_alloc(lx, ux);
120 for (i = ly; i < uy; ++i)
121 {
122 for (j = lx; j < ux; ++j)
123 {
124 for (k = imptrlz; k < imptruz; k++)
125 {
126 row = imptrs[k][i];
127 data[k] = row[j];
128 }
129
130 y = im_get_pixf(im_slope, i, j);
131
132 if (select1 == select2)
133 {
134 a = 1.0;
135 d = 0.0;
136 tissue = m->normal[select1]*pure_density(m,&data[imptrlz], select1);
137 }
138 else
139 {
140 a = part_fraction(m,&data[imptrlz],select1,select2);
141 fract1 = part_density(m,&data[imptrlz],select1,select2,a);
142 fract2 = part_density(m,&data[imptrlz],select2,select1,(1.0-a));
143 if (fract1+fract2 > 0.0) a = fract1/(fract1+fract2);
144 d = mD_dist_calc(m,select1,select2);
145 tissue = m->par_dens[select1][select2]*fract1;
146 }
147
148
149 Pgrad =mD_grad_func(a,d,y,select1,select2,Sk,s0,&sf);
150 row1[j] = tissue*Pgrad;
151
152 }
153 im_put_rowf(row1, im_prob, i, lx, ux);
154 }
155 fvector_free(row1,lx);
156 fvector_free(data,imptrlz);
157
158 return (im_prob);
159 }
160
161
162 Imrect *prob_im( Mixmodel *m, Imrect *im_prob, Imregion *roi, int select1, int select2)
163 {
164 int lx, ux, ly, uy;
165 int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
166 double tissue,a;
167 float *row, *row1, *data;
168 int i,j,k;
169 void ***imptrs = NULL;
170
171 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly,
172 &imptruy, &imptrlx, &imptrux);
173 data = (float *)fvector_alloc(imptrlz, imptruz);
174
175 if (im_prob == NULL)
176 {
177 error("Density image memory not allocated!", warning);
178 return(NULL);
179 }
180
181 if (m == NULL)
182 {
183 error("Input file is missing!", warning);
184 return(NULL);
185 }
186 if (imptrlx > (lx = roi->lx)) lx = imptrlx;
187 if (imptrux < (ux = roi->ux)) ux = imptrux;
188 if (imptrly > (ly = roi->ly)) ly = imptrly;
189 if (imptruy < (uy = roi->uy)) uy = imptruy;
190
191 row1 = fvector_alloc(lx, ux);
192 for (i = ly; i < uy; ++i)
193 {
194 for (j = lx; j < ux; ++j)
195 {
196 for (k = imptrlz; k < imptruz; k++)
197 {
198 row = imptrs[k][i];
199 data[k] = row[j];
200 }
201
202 if (select1 == select2)
203 tissue = m->normal[select1]*pure_density(m,&data[imptrlz], select1);
204 else
205 {
206 a = part_fraction(m,&data[imptrlz],select1,select2);
207 tissue = m->par_dens[select1][select2]*
208 part_density(m,&data[imptrlz],select1,select2,a);
209 }
210 row1[j] = tissue;
211
212 }
213 im_put_rowf(row1, im_prob, i, lx, ux);
214 }
215 fvector_free(row1,lx);
216 fvector_free(data,imptrlz);
217 return (im_prob);
218 }
219
220
221 Imrect *prob_im_tot(Mixmodel *m, void ***pim_arr,Imregion *roi)
222 {
223 Imrect *im,*im_tot1, *im_tot2,*im_bkgd;
224 int i,j;
225 float n, current;
226 void ***imptrs = NULL;
227 int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
228 int lx, ux, ly, uy;
229
230 if (pim_arr == NULL)
231 {
232 error("Probability densities not calculated!", warning);
233 return(NULL);
234 }
235 if (m == NULL)
236 {
237 error("Input file is missing!", warning);
238 return(NULL);
239 }
240
241 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
242 &imptrux);
243
244 if (imptrlx > (lx = roi->lx)) lx = imptrlx;
245 if (imptrux < (ux = roi->ux)) ux = imptrux;
246 if (imptrly > (ly = roi->ly)) ly = imptrly;
247 if (imptruy < (uy = roi->uy)) uy = imptruy;
248
249 im_tot1 = im_alloc(imptruy-imptrly, imptrux-imptrlx, roi, float_v);
250 im_bkgd = im_alloc(imptruy-imptrly, imptrux-imptrlx, roi, float_v);
251
252 for(i = 0; i < m->nmix; i++)
253 {
254 for(j = 0; j < m->nmix; j++)
255 {
256 if (pim_arr[i][j] != NULL)
257 {
258 /*get the memory address for the image*/
259 im = (Imrect *)pim_arr[i][j];
260 im_tot2 = im_sum(im_tot1,im);
261 im_free(im_tot1);
262 im_tot1=im_tot2;
263 }
264 }
265 }
266 n = (ux-lx) *(uy-ly);
267 for (i = ly; i < uy; i++)
268 {
269 for (j = lx; j < ux; j++)
270 {
271 current = im_get_pixf(im_tot1,i,j);
272 if (current< n*m->background)
273 im_put_pixf(n*m->background-current,im_bkgd,i,j);
274 }
275 }
276 im_tot2 = im_sum(im_tot1,im_bkgd);
277 im_free(im_tot1);
278 im_free(im_bkgd);
279 return (im_tot2);
280 }
281
282 Imrect *prob_im_norm(Imrect *im_orig, Imrect *im_tot)
283 {
284 /* Normalisation of each segmented image */
285 Imrect *im_norm;
286 float threshold = 0.0;
287 float rep_val = 0.0;
288
289 im_norm = im_div(im_orig,im_tot,threshold,rep_val);
290 return (im_norm);
291 }
292
293
294 /*
295 * Commented out - GAB 19 Feb 2003
296 * Not completely deleted because it may be useful for testing.
297 *
298 List *prob_all(Mixmodel *m, double **select_matrix, void *** pim_array)
299 {
300 Imrect *im_tot1,*im_tot2,*im, *im_tot;
301 Imrect *im_seq;
302 int imptrlx, imptrux, imptrly, imptruy;
303 int loop,loop1;
304 List *lptr = get_current_seq_current_el();
305 Imregion *roi;
306
307 if ((im_seq = (Imrect *)(lptr->to))== NULL)
308 {
309 error("Image sequence is missing!", warning);
310 return;
311 }
312 roi = im_seq->region;
313 if (roi == NULL)
314 return;
315 if (m == NULL)
316 {
317 error("Input file is missing!", warning);
318 return;
319 }
320 if (pim_array == NULL)
321 {
322 error("Probability densities not calculated!", warning);
323 return;
324 }
325
326 im_tot = prob_im_tot(m,pim_array,roi);
327
328 for(loop = 0; loop < m->nmix; loop++)
329 {
330 im_tot1 = im_alloc(imptruy-imptrly, imptrux-imptrlx, roi, float_v);
331 for(loop1 = 0; loop1 < m->nmix; loop1++)
332 {
333 if (pim_array[loop][loop1] != NULL)
334 {
335 im = (Imrect *)pim_array[loop][loop1];
336 im = prob_im(m,im,roi,loop,loop1);
337 im_tot2 = im_sum(im_tot1,im);
338 im_free(im_tot1);
339 im_tot1=im_tot2;
340 }
341 }
342 im_tot2 = prob_im_norm(im_tot1,im_tot);
343 stack_push(im_copy(im_tot2), IMRECT, im_free);
344 im_free(im_tot1);
345 }
346 im_free(im_tot);
347 im_free(im_tot2);
348
349 }
350 *
351 * End of "commented out" for function prob_all - GAB 19 Feb 2003
352 */
353
354 Imrect *prob_class(Mixmodel *m, void *** pim_array, char *str_name)
355 {
356 Imrect *im_tot1=NULL, *im_tot2=NULL, *im=NULL, *im_tot=NULL;
357 Imrect *im_seq=NULL;
358 void ***imptrs = NULL;
359 int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
360 int loop,loop1,i;
361 List *lptr = get_current_seq_current_el();
362 Imregion *roi=NULL;
363 int str_check=0;
364
365 if ((im_seq = (Imrect *)(lptr->to))== NULL)
366 {
367 error("Image sequence is missing!", warning);
368 return(NULL);
369 }
370 roi = im_seq->region;
371
372 if (roi == NULL)
373 return(NULL);
374 if (m == NULL)
375 {
376 error("Input file is missing!", warning);
377 return(NULL);
378 }
379
380 if (pim_array == NULL)
381 {
382 error("Probability densities not calculated!", warning);
383 return(NULL);
384 }
385 if(str_name == NULL)
386 {
387 error("Tissue type not specified \n", warning);
388 return(NULL);
389 }
390
391 for (i=0;i<m->nmix;i++)
392 {
393 if(strcmp(m->tissue_type[i],str_name) == 0) str_check =1;
394 }
395 if( str_check == 0)
396 {
397 error("Wrong tissue type specified! \n", warning);
398 return(NULL);
399 }
400
401 im_tot = prob_im_tot(m,pim_array,roi);
402
403 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
404 &imptrux);
405
406 for(loop = 0; loop < m->nmix; loop++)
407 {
408 if( strcmp(str_name, m->tissue_type[loop]) == 0)
409 {
410 im_tot1 = im_alloc(imptruy-imptrly, imptrux-imptrlx, roi, float_v);
411 for(loop1 = 0; loop1 < m->nmix; loop1++)
412 {
413 if (pim_array[loop][loop1] != NULL)
414 {
415 im = (Imrect *)pim_array[loop][loop1];
416 im_tot2 = im_sum(im_tot1,im);
417 im_free(im_tot1);
418 im_tot1=im_tot2;
419 }
420 }
421 im_tot2 = prob_im_norm(im_tot1,im_tot);
422 im_free(im_tot1);
423 }
424 }
425 im_free(im_tot);
426 return(im_tot2);
427 }
428
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.