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

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

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