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

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

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

  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 

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